Introduction

Who is the instructor?

Materials

## Required packages in this workshop
lib2install <- c("metaSEM", "lavaan", "semPlot", "readxl")

## Install them automatically if they are not available on your computer
for (i in lib2install) {
  if (!(i %in% rownames(installed.packages()))) install.packages(i)
}

Objectives of the workshop

  • This workshop provides a practical introduction of meta-analytic structural equation modeling (MASEM) using the open-source R statistical platform.
  • Expectations:
    1. Learn essential ideas of SEM and meta-analysis;
    2. Learn basic and some advanced techniques in MASEM;
    3. Learn how to conduct and interpret MASEM in R.
    4. But don’t expect to be an expert in MASEM in 3 hours!

A crash course in SEM

What is SEM?

  • SEM represents a family of related techniques.
  • It is also known as covariance structural analysis, covariance structure model, analysis of covariance structures, analysis of correlation structure, LISREL model (in the old days), etc.

Key concept 1: Observed vs. latent variables

  • Latent variables:
    • Abstract and hypothetical constructs.
    • They cannot be measured directly.
    • Most constructs in psychology and social sciences are unobservable.
    • For example, motivation, stress, depression, intelligence, and satisfaction.
  • Observed or measured variables:
    • Indicators of the latent constructs.
    • For example, items to measure depression, test scores of intelligence.

An example

  • Satisfaction with life scale (SWLS)
    1. In most ways my life is close to my ideal.
    2. The conditions of my life are excellent.
    3. I am satisfied with my life.
    4. So far I have gotten the important things I want in life.
    5. If I could live my life over, I would change almost nothing.
  • Physical health, e.g.,
    1. Blood pressure
    2. Blood sugar level
    3. Cholesterol levels

Key concept 2: Models of the relationship among the constructs

  • Most statistical techniques, e.g., regression, are limited to one dependent variable (DV).
  • SEM allows researchers to test models with a complicated relationship.
    • Models are usually represented by path diagrams.
  • Model fitting:
    • Overall model fit: Does the proposed model fit the data?
    • Individual fit: Is the parameter estimate statistically significant?

Key concept 3: Notations for path diagrams

  • Rectangles or squares: observed (or measured) variables
  • Circles or ellipses: unobserved (or latent) variables
  • Triangles: means or intercepts (may not be used if we are not interested in the mean structure)
  • Direct arrows: predictions or “causes”
  • Double arrows between two variables: covariances
  • Double arrows on the same variables: variances

Path analysis

  • Linear relationships among observed variables.
  • No latent variable.
  • Multiple regression may also be used to fit this model.
  • For example, we can estimate an indirect effect by constructing a function of the parameters (\(a*b\)). We may compute the 95% confidence interval (CI) on the indirect effect.

Confirmatory factor analysis (CFA)

  • Linear relationships among latent and observed variables.
  • No direct effect on the latent variables.
  • It is used to test the construct validity of psychological constructs.

SEM

  • It combines both CFA and path analysis.
  • There are observed and latent variables.
  • It may include direct effects among the latent variables.

How to specify a structural equation model?

Model specification

  • We need to instruct the computer programs on what models we would like to fit.
  • Approaches to specify SEM
    • Matrices: LISREL and OpenMx (open0source software)
    • Equations: Mplus, EQS, and lavaan (open-source software)
    • Path diagrams: Amos

Model-identification (1)

  • After specifying a model, we need to evaluate whether the proposed model is testable.
  • Model identification concerns whether there is a unique solution for the model being tested.
  • If the model is not identified, there is no solution for the proposed model. That is, we are asking a question that has no answer.
  • Degrees of freedom (dfs) of a model:
    • No. of pieces of information available on the covariance matrix: \(p^*=p(p+1)/2\) where \(p\) is the no. of variables.
    • No. of free parameters in the model: \(q\)
    • \(df = p^*-q\)

Model-identification (2)

  • A necessary but not sufficient condition for identifying any SEM model is df >= 0.
  • Underidentification:
    • df < 0
    • no solution
    • e.g., 3=x+y
  • Just identification:
    • df =0
    • perfect fit
    • e.g., 3=x+y and 1=x-y then x=2 and y=1
  • Overidentification:
    • df > 0
    • no perfect solution
    • e.g., 3=x+y, 2=x-y and 5=2x+y

Model-identification (3)

  • All regression models are just identified models (df=0).
  • We cannot empirically test whether or not the proposed regression models are correct. (They are always correct by definition.)
  • Over-identified models (dfs > 0) are usually preferred in CFA/SEM.
  • Using a falsification approach, we compare the proposed model (theory) against the data.

Data in the analysis

  • If the data are multivariate normal, we may use the covariance matrix as the input for data analysis.
  • The above assumption may simplify the computations.
  • Degrees of freedom (dfs) of the models can be calculated and checked easily.
  • When there are missing data or the data are not continuous, raw data are required.
  • Using the simple regression as an example:
    • \(S = \left[ {\begin{array}{cc} Var(y) & \\ Cov(y,x_1) & Var(x_1) \\ \end{array} } \right] = \left[ {\begin{array}{cc} 10.61 & \\ 5.54 & 7.78\\ \end{array} } \right]\)
  • The no. of pieces of information available: \(p^* = p(p+1)/2=3\), where p is the no. of observed variables.
  • But we only have the summary statistics, e.g., correlation matrices, in MASEM.
Sample data
x1 y
7.49 5.27
10.66 5.29
9.61 4.40
14.43 7.22
10.58 5.37
11.59 4.91
7.09 3.83
13.57 6.53
5.87 2.26
8.20 11.93

Model specification

  • We can specify a simple regression via a path diagram. We skip the means for the time being.

  • The model in the equation: \(y = b_0 + b_1 x_1 + e_y\)
  • Given this model, we can calculate the model implied covariance matrix:
    • \(\Sigma(\theta) = \left[ {\begin{array}{cc} b^2_1 Var(x_1) + Var(e_y) & \\ b_1 Var(x_1) & Var(x_1) \\ \end{array} } \right]\)
  • The above equation (in plain language) means that the population covariance matrix \(\Sigma\) is a function of \(\theta\).

Parameter estimation (1)

  • Under the proposed model, the model implied covariance matrix is \(\Sigma=\Lambda \Phi \Lambda^T+\Psi\).
  • We try to find the parameter estimates such that the model implied covariance matrix (theory) \(\hat{\Sigma}\) is as close to the sample covariance matrix (data) \(S\) as possible, i.e., \(\hat{\Sigma} \approx S\).
  • We usually use the maximum likelihood (ML) estimation method.

Parameter estimation (2)

  • Data = Model + Error;
  • Date: Sample covariance matrix (S) vs. Model: Model implied covariance matrix (\(\hat{\Sigma}\)).
  • For example, \(S = \left[ {\begin{array}{c|cccc} & x1 & x2 & x3 & x4\\ \hline x1& 4& & & \\ x2& 2.0& 3& & \\ x3& 2.3& 3.4& 5 &\\ x4& 2.4& 3.2& 3.0 &6\\ \end{array} } \right]\) and \(\hat{\Sigma} = \left[ {\begin{array}{c|cccc} & x1 & x2 & x3 & x4\\ \hline x1& 3.4& & & \\ x2& 1.5& 3.2& & \\ x3& 2.0& 3.2& 5.2 &\\ x4& 2.6& 3.5& 3.2 &5.6\\ \end{array} } \right]\).
  • Error = Data - Model = \(S-\hat{\Sigma}\).
  • The amount of Error may indicate how good the model fits the data.

Model testing

Model evaluation in SEM

  • There are two major tasks in model evaluation.
  • Overall model fit: testing whether the proposed model as a whole fits the data.
  • Individual parameter estimates: testing whether the parameter estimates are significant.
  • Note. If the overall model does not fit the data, we do not test and interpret the parameter estimates.

Chi-square test statistics

  • Chi-square test (also known as the likelihood ratio (LR) test):
    • If the proposed model is correct, the test statistic has a chi-square distribution with (p*-q) df.
    • This is a “badness-of-fit” index: large chi-square statistic indicates a poor fit.
    • The proposed model is rejected at .05 if the test statistic is larger than the critical value.

Issues of chi-square test statistics

  • Violation of underlying assumptions:
    • Data (or residues) are normally distributed.
    • Large samples are required.
    • When data are not normally distributed, especially in clinical studies, or small sample sizes (e.g., N=100 or 200), the test statistic may not follow a chi-square distribution.
  • Sensitive to sample size:
    • All proposed models will be rejected when the sample sizes are large enough.
    • Large samples work against researchers!

Goodness-of-fit indices

  • Many SEM users are aware of the problems associated with the chi-square test.
  • There are many goodness-of-fit indices developed as alternative measures.
  • There are many goodness-of-fit indices in the market!

Incremental fit indices

  • They measure the relative improvement in fit by comparing the target model against the baseline model.
  • The target model is usually the proposed model.
  • The baseline model is usually the model stating that all variables are uncorrelated. It is known as the independence model. It can be considered the worst model.
  • Normed fit index (NFI):
    • It measures the proportionate reduction in the chi-square values from the baseline model to the hypothesized model.
  • Non-normed fit index (NNFI) which is known as Tucker-Lewis index (TLI):
    • It is similar to NFI but adjusts the complexity of the model.
  • Comparative fit index (CFI):
    • Makes sure that the maximum is 1.
  • What is a well-fitted model?
    • Conventional rule of thumb (without any empirical support): at least > 0.9

Residual-based indices (1)

  • When the model fits well, the residuals (the difference between the model implied covariance matrix and the sample covariance matrix) should be small.
  • Standardized root mean square residual (SRMR)
    • It measures the average value of the standardized residuals.
    • It ranges from zero (perfect fit) to one (very poor fit).
    • Rule of thumb: A well-fitted model < .05.

Residual-based indices (2)

  • Root mean square error of approximation (RMSEA) (Browne & Cudeck, 1993)1
    • Similar to SRMR.
    • Advantage: Confidence intervals on RMSEA are available on most SEM packages.
    • Rules of thumb:
      • Close fit: < 0.05
      • Reasonable fit: 0.05 - 0.08
      • Inadequate fit: > 0.1

What do we need to report in research articles?

  • We usually report the chi-square test statistic, and it associated df and p-value, some incremental fit indices, and some residual-based indices.
  • What is a well-fitted model? One popular approach is the combinational rules (Hu & Bentler, 1999):2
    • NNFI (TLI) or CFI > 0.95 and SRMR < .09 OR RMSEA < .05 and SRMR < .06
  • Although this recommendation has been widely applied, it is not without criticisms (Marsh, Hau, & Wen, 2004).3

A crash course in meta-analysis

Are systematic review and meta-analysis the same?

  • Systematic review:
    • It aims to provide a comprehensive literature search with pre-defined eligibility criteria.
    • It focuses on minimizing bias in a literature review by using pre-defined criteria to replicate the literature search.
    • The main question is how to identify the relevant studies correctly.
  • Meta-analysis:
    • It combines and synthesizes findings with statistical models.
    • The main question is how to analyze the quantitative data properly.

Not all types of evidence are equal.

  • There is a hierarchy of evidence in medical science (e.g., Guyatt et al., 1995),4
    • Experimental study: randomized control trials (participants are randomly assigned to vaccine vs. non-vaccine groups)
    • Observational study:
      • Cohort studies: recruit participants with and without vaccines and follow them to see if they will get COVID-19 in the future.
      • Case-control studies: recruit participants with and without COVID-19 and trace whether they have been vaccinated.
  • In social sciences, we also have similar ideas, e.g.,
    • Meta-analysis > experimental studies > observational studies > case studies.

Problems in empirical research

  • An artificial example of the correlations between income and self-esteem:

    Study Sample size (ni) Correlation (ri) p value
    1 26 .13 .53
    2 42 .35 .03
    3 20 -.10 .68
    4 40 .30 .07

  • If you were the researchers conducting Studies 2 and 4, you might argue that there was a positive relationship.

  • If you were the researchers conducting Studies 1 and 3, you might conclude that there was no linear relationship.

  • The results look contradictory to each other.

Problems of individual study

  • Measurement artifacts, e.g., the unreliability of measurement, restriction of range → significant or non-significant results?
  • Small sample sizes → significant or non-significant results?
  • Low statistical power → significant or non-significant results?
  • Questionable research practices, e.g., reporting significant results, multiple testing without controlling for it → significant or non-significant results?
  • Questions:
    • Can we still trust the published findings if there are so many issues in the published results?
    • If not, what are the alternatives?

Objectives of a meta-analysis

  • Draw general conclusions on a particular topic.
  • Test the homogeneity (consistency) of the existing findings.
  • Account for the heterogeneity of effect sizes.
  • Estimate an average effect size.
  • Test potential moderators if the studies are heterogeneous.
  • Assess robustness of the findings, e.g., publication bias.

Effect sizes and their sampling variances

What is an effect size? (1)

  • Effect sizes and their sampling variances are the building blocks of meta-analysis.
  • There are 3 considerations:
    1. Strength
    2. Direction
    3. Independence of sample size
  • Questions:
    1. Can a t statistic be used as an effect size in meta-analysis?
    2. Can an \(R^2\) and \(\Delta R^2\) in multiple regression be used as effect sizes in meta-analysis?
    3. Can an \(\eta^2\) and \(\omega^2\) in ANOVA be used as effect sizes in meta-analysis?
    4. Can a regression coefficient \(\hat{\beta}\) be used as an effect size in meta-analysis?

What is an effect size? (2)

  • We use y and v to represent a generic effect size and its sampling variance. y can be a correlation coefficient, a standardized mean difference, or other effect sizes.
  • Once we have calculated the effect sizes and their sampling variances, we may apply the same meta-analytic procedures to the data.

Common families of effect sizes

  • There are three popular families:
    • r: Pearson correlation, which is popular in observational studies.
    • d: (standardized) mean difference, which is popular in experimental designs (control vs. intervention groups).
    • OR: odds ratio, which is popular in medical research with binary outcomes.
Cheung and Vijayakumar (2016, Table 1)
Cheung and Vijayakumar (2016, Table 1)

Correlation (r)

  • Raw correlation
    • Correction on artifacts may be applied (Schmidt & Hunter, 2015).5
    • \(y_r = r\)
    • \(v_r = SE^2_{r} =\frac{(1-r^2)^2}{n-1}\)
  • Fisher’s z transformed score
    • It may stabilize the sampling distribution of the correlation coefficient (Hedges & Olkin, 1985)6.
    • \(y_z = 0.5*\mathrm{log}\frac{1+r}{1-r}\)
    • \(v_z = SE^2_{z} = \frac{1}{n-3}\)
    • One nice feature of the Fisher’s z score is that its sampling variance is independent of z. Thus, any sampling error in estimating r (or z) does not affect its estimated sampling variance.
    • After the meta-analysis, the average \(y_z\) is converted back to \(y_r=\frac{e^{2z}+1}{e^{2z}-1}\) for ease of interpretations.

Models for meta-analysis

  • Generally, there are two models for meta-analysis—fixed- and random-effects models (e.g., Hedges & Vevea, 1998).7
  • Fixed-effect or common-effect model
    • The population effect sizes are assumed to be equal for all studies
    • The differences in the observed effect sizes are due to sampling error.
    • We can generalize the findings similar to those included in the meta-analysis regarding research design, measures, and samples.
    • We cannot generalize the findings beyond the studies included in the meta-analysis.
  • Random-effects model
    • The population effect sizes can be different across studies.
    • The differences in the observed effect sizes are due to a combination of true difference (variance component) and sampling error.
    • We can generalize the findings beyond the studies included in the meta-analysis.

Random-effects model

  • The population effect sizes can differ across studies (Hedges & Vevea, 1998).8
  • The differences in the observed effect sizes are due to a combination of true difference (variance component) and sampling error.
  • We can generalize the findings beyond the studies included in the meta-analysis.
  • It is usually preferred in meta-analysis.
  • Model: \(y_i = \beta_\mathrm{R} + u_i + e_i\),
    • where \(y_i\) is the observed effect size, \(\beta_\mathrm{R}\) is the average population effect size under a random-effects model, and \(\mathrm{Var}(e_i)=v_i\) is the conditional sampling variance of \(y_i\), and \(u_i\) is the study-specific random effect.
    • \(\tau^2=\mathrm{Var}(u_i)\) is the variance of the random effects. It indicates the between-study heterogeneity.
    • The fixed-effect model is a special case of the random-effects model when \(\tau^2=0\).

Differences between a fixed-effect and random-effects model

  • Borenstein et al. (2009)9 have provided illustrations between fixed-effect and random-effects models.

Fixed-effect model

  • All studies share the same common population effect \(\beta_\mathrm{F}\) (circles).
Borenstein et al. (2009, Figure 11.1)
Borenstein et al. (2009, Figure 11.1)
  • The observed difference in the sample effect sizes \(y_i\) (squares) are only due to the sampling error with a variance of \(v_i\).
Borenstein et al. (2009, Figure 11.3)
Borenstein et al. (2009, Figure 11.3)

Random-effects model

  • There is an average population effect \(\beta_\mathrm{R}\) (triangle) with a distribution (Var=\(\tau^2\))
Borenstein et al. (2009, Figure 12.1)
Borenstein et al. (2009, Figure 12.1)
  • Each study has its own true population effect \(\beta_i\) (circles), which are not observed by us. That is, the true population effect is latent.
Borenstein et al. (2009, Figure 12.2)
Borenstein et al. (2009, Figure 12.2)
  • The sample effect sizes \(y_i\) (squares) are conditionally sampled with \(v_i\) due to the sampling error.
Borenstein et al. (2009, Figure 12.4)
Borenstein et al. (2009, Figure 12.4)

Issues in interpreting \(\tau^2\) (1)

  • The value of \(\tau^2\) cannot be compared across effect sizes because it depends on the type of effect size.
  • Higgins and Thompson (2002)10 suggested three criteria for measuring heterogeneity:
    • Dependence on the extent of heterogeneity
    • Scale invariance (\(\tau^2_r\) and \(\tau^2_d\) should be comparable)
    • Size invariance (k; it does not depend on the number of studies)
    • Note. It does not say that it is sample size (N) invariance.

Issues in interpreting \(\tau^2\) (2)

  • They further proposed several indices to quantify the degree of heterogeneity,
    • \(I^2=\frac{Q-(k-1)}{Q}\).
    • When the estimated \(I^2\) is negative, it is truncated to 0.
  • An alternative way to express it is,
    • \(I^2=\frac{\hat{\tau}^2}{\hat{\tau}^2 + \tilde{v_i}}\), where \(\tilde{v_i}\) is the “average” sampling variance.
  • It can be interpreted as the proportion of total variation due to heterogeneity between studies.
  • As a rule of thumb, \(I^2\) of 25%, 50%, and 75% can be considered as low, moderate, and high heterogeneity. This rule of thumb was based on meta-analyses of medical research. They may or may not be the same in other disciplines.
  • When the sample sizes in the primary studies are getting larger and larger, \(I^2\) approaches 1.
  • Note. \(\tau^2\) is an absolute measure of heterogeneity, whereas \(I^2\) is a relative measure of heterogeneity depending on \(\tilde{v_i}\) (Borenstein et al., 2017).11 It is advisable to report both indices.

Confidence interval

  • Confidence interval (CI):
    • A 95% confidence level represents the theoretical long-run frequency (i.e., 95% proportion) of confidence intervals that contain the true value of the unknown population parameter.
    • CIs can only be formed on parameters, e.g., \(\beta_\mathrm{F}\), \(\beta_\mathrm{R}\), and \(\tau^2\).
    • For example, if the \(\hat{\beta}_\mathrm{R}=.2\) and 95% CI = (.1, .3), \(\hat{\beta}_\mathrm{R}\) is statistically significant at \(\alpha=.05\).

Why are researchers concerned about the heterogeneity of studies?

  • Do you prefer homogeneous or heterogeneous findings? Why?
  • What are the implications of homogeneous findings?
  • What are the implications of heterogeneous findings?

Exploring heterogeneity

  • Why the effect sizes are heterogeneous?
    • Sampling error (bad luck!)
    • Measurement artifacts, e.g., unreliability, range restriction, differences in samples, measures, etc.
    • Presence of moderators, e.g., samples, measures, etc.
  • When there is excess heterogeneity, we may want to model the effect sizes with the study characteristics.
  • A subgroup analysis and mixed-effects model are the common technique to do it.

Subgroup analysis

  • If the moderator is a categorical variable, e.g., experimental vs. observational, we may run a separate meta-analysis on each level.
  • Advantages:
    • It is simple.
    • Different groups may have different estimates of heterogeneity variances.
  • Disadvantages:
    • It is slightly more complicated if we want to test whether the categorical variable is statistically significant.
    • It does not provide an \(R^2\) on the moderator.
  • Readers may refer to Rubio-Aparicio et al. (2020)12 for a comparison between the subgroup analysis and meta-regression in handling categorical moderators.

Mixed-effects model

  • Mixed-effects model or meta-regression may be used to explore the heterogeneity (Konstantopoulos & Hedges, 2004).13
  • The study characteristics are known as the covariates, predictors, or moderators in meta-analysis.

What is MASEM?

MASEM

  • Meta-analysis may also be used to synthesize findings in structural equation models. The class of techniques is called meta-analytic structural equation modeling (MASEM) (Cheung, 2021).14
  • It uses meta-analytic techniques to combine correlation matrices and fits structural equation models on the average correlation matrix. It has the best of both worlds.

Benefits of MASEM to SEM

  • MASEM allows researchers to integrate findings over several studies.
  • MASEM may test models without testing in individual studies.
  • MASEM may study how the model varies according to the study characteristics via moderator testing.

Benefits of MASEM to meta-analysis

  • MASEM tests theoretically relevant models, e.g., a mediation model, rather than the individual effects.
  • MASEM allows researchers to compare several models.

Problems in primary research using SEM

  • Low statistical power in psychological studies;
  • Contradictory findings with significant tests;
  • Confirmation bias: reluctant to consider alternative models in SEM;
  • Different models have been proposed by various researchers.
  • Conducting more empirical studies may not necessarily decrease the uncertainty of a particular topic if the findings are inconsistent (National Research Council, 1992).15
  • MASEM may be used to address many of these issues.

Terms used in the literature

  • Several names have been used interchangeably in the literature,
    • meta-analytic path analysis;
    • meta-analysis of factor analysis;
    • meta-analytical structural equations analysis;
    • path analysis of meta-analytically derived correlation matrices;
    • SEM of a meta-analytic correlation matrix;
    • path analysis based on meta-analytic findings; and
    • model-based meta-analysis.
  • We use the generic term meta-analytic structural equation modeling (MASEM) to describe this class of techniques.

Example: Brown and Stayman (1992)

  • Topic: Antecedents and consequences of attitudes toward the advertisement (Ad)16
    • No. of variables: 5
    • No. of studies: 47
    • Pooled sample size across studies: +4,600
Brown and Stayman (1992)
Brown and Stayman (1992)

Example: Premack and Hunter (1988)

  • Topic: Individual unionization decisions17
    • No. of variables: 6
    • No. of studies: 14
    • Pooled sample size across studies: +2,800
Premack and Hunter (1988)
Premack and Hunter (1988)

Example: Norton et al. (2013)

  • There are ten different factor structures on the Hospital Anxiety and Depression Scale supported by some empirical data.18
  • The authors identified 28 independent samples from 21 studies (N=21,820).
  • They found that the bifactor structure consisting of a general distress factor and anxiety and depression group factors fitted the data best.
Norton et al. (2013)
Norton et al. (2013)

Approaches to MASEM

  • Nearly all methods use a two-stage approach to conducting MASEM:
    • Stage 1 analysis: Combine the correlation matrices into a pooled correlation matrix;
    • Stage 2 analysis: Fit structural equation models on the pooled correlation matrix.
  • Common approaches:
    • Univariate approach (Hunter & Schmidt, 1990; Viswesvaran & Ones 1995);19
    • Generalized least squares (GLS; Becker, 1992);20
    • Two-stage structural equation modeling (TSSEM; Cheung, 2014; Cheung & Chan, 2005).21, 22
    • One-stage MASEM (OSMASEM; Jak & Cheung, 2020).23
Jak and Cheung (2020, Table 1)
Jak and Cheung (2020, Table 1)

Univariate approach (Viswesvaran & Ones, 1995)

  • Stage 1 analysis
    • It synthesizes the correlation coefficients in a correlation matrix as if the correlation coefficients were independent.
    • Incomplete correlation coefficients are handled by either listwise deletion or pairwise deletion (more popular).
    • Problem: when the pairwise deletion is used, missing correlations are assumed to be missing completely at random (MCAR).
  • Stage 2 analysis
    • The pooled correlation matrix is used as if it were a covariance matrix.
    • Researchers usually use the harmonic mean of the sample sizes as the sample size.

Problems of the univariate approach

  • The correlation matrix is analyzed as if it were a covariance matrix. The test statistics and the SEs may be incorrect (Cudeck, 1989).24
  • One sample size is used to represent the precision of the average correlation matrix.
  • The univariate approaches are based on the assumption that MASEM=MA+SEM.
  • For example, \(R_1 = \left[ {\begin{array}{c|ccc} & x1 & x2 & x3\\ \hline x1& 1.0& & & \\ x2& 0.6& 1.0& & \\ x3& 0.5& 0.7& 1.0 &\\ \end{array} } \right]\), \(R_2 = \left[ {\begin{array}{c|ccc} & x1 & x2 & x3\\ \hline x1& 1.0& & & \\ x2& 0.6& 1.0& & \\ x3& NA& NA& NA&\\ \end{array} } \right]\), \(R_3 = \left[ {\begin{array}{c|ccc} & x1 & x2 & x3\\ \hline x1& 1.0& & & \\ x2& NA& 1.0& & \\ x3& 0.7& NA& 1.0 &\\ \end{array} } \right]\), and \(R_4 = \left[ {\begin{array}{c|ccc} & x1 & x2 & x3\\ \hline x1& 1.0& & & \\ x2& NA& 1.0& & \\ x3& NA& 0.5& 1.0 &\\ \end{array} } \right]\).
  • The average correlation matrix is \(\bar{R} = \left[ {\begin{array}{c|ccc} & x1 & x2 & x3\\ \hline x1& 1.0& & & \\ x2& 0.6& 1.0& & \\ x3& 0.6& 0.6& 1.0 &\\ \end{array} } \right]\).

One (sample) size does not fit all!

  • There are some problems fitting SEMs with a single sample size in the univariate approach.
  • Let us consider the average correlation matrix as an example.
  • The sample sizes for the correlations \(r_{21}\), \(r_{31}\), and \(r_{32}\) are 200, 500, and 1,000, respectively.
  • The arithmetic mean, the harmonic mean, and the median are 567, 375, and 500, respectively.
  • If the harmonic mean (375) is used as the sample size in SEM, some estimated SEs are over-estimated while others are under-estimated.
  • The sample size affects the chi-square test, some goodness-of-fit indices, and SEs.
  • \(\bar{R} = \left[ {\begin{array}{c|ccc} & x1 & x2 & x3\\ \hline x1& 1.0& & & \\ x2& 0.6& 1.0& & \\ x3& 0.6& 0.6& 1.0 &\\ \end{array} } \right]\).

Problems of treating a correlation matrix as a covariance matrix

  • Analysis of the correlation matrix:
    • The diagonals are always 1. They do not carry any information.
    • The diagonals of \(\hat{\Sigma}\) are precisely 1.
  • The chi-square test or SEs can be incorrect if we treat a correlation matrix as if it was a covariance matrix in SEM.
  • Cudeck (1989) warned of the problems of treating correlation matrices as covariance matrices in SEM almost 30 years ago. However, many researchers still ignore the warnings in MASEM today.25

A simulation study comparing some of these methods

  • Jak and Cheung (2020)26 (see Jak & Cheung (2022))27 conducted a simulation study comparing several approaches.
    • The parameter estimates were similar among all approaches.
    • GLS, TSSEM, and OSMASEM had similar performance, whereas the univariate approach did not perform well.
Jak and Cheung (2020)
Jak and Cheung (2020)
Jak and Cheung (2020)
Jak and Cheung (2020)

TSSEM approach

  • The TSSEM approach is very similar to the GLS approach, but it is based on SEM.
  • Fixed-effects TSSEM (Cheung & Chan, 2005)28
    • Researchers are only interested in the studies included in the meta-analysis;
    • A common correlation/covariance matrix is assumed.
  • Random-effects TSSEM (Cheung, 2014)29
    • Studies are samples of a larger population;
    • Researchers may want to generalize the findings beyond the studies included;
    • Studies may have their own correlation/covariance matrices.

Random-effects TSSEM

  • The fixed-effects TSSEM has been extended to a random-effects model (Cheung, 2014), which is mathematically the same as the random-effects GLS proposed by Becker (1992).
  • There are several pros and cons of applying a random-effects model:
    • Pros: Treating differences across studies as random effects;
    • Cons: Difficult to explain study differences by using study characteristics as the moderators.
  • Stage 1 analysis:
    • Multivariate random-effects model is used to synthesize correlation matrices;
  • Stage 2 analysis:
    • WLS estimation is used to fit structural equation models on the average correlation matrix.

A conceptual overview of TSSEM

Jak and Cheung (2020, Figure 1)
Jak and Cheung (2020, Figure 1)

Random-effects TSSEM: Stage 1 analysis (1)

  • A multivariate meta-analysis is used to synthesize the correlation matrices
    • \(\boldsymbol{r}_i = \boldsymbol{\rho}_\mathrm{R} + \boldsymbol{u}_i + \boldsymbol{e}_i\),
    • where \(\boldsymbol{e}_i \sim \mathcal{N}(\boldsymbol{0}, V_i)\) is the known sampling covariance matrix and \(\boldsymbol{u}_i \sim \mathcal{N}(0, T^2)\) is the heterogeneity covariance matrix that has to be estimated.
  • The SEM-based multivariate meta-analysis is used to obtain the average correlation matrix \(R_\mathrm{R}\) and its asymptotic sampling covariance matrix \(V_\mathrm{R}\).

Random-effects TSSEM: Stage 1 analysis (2)

  • \(I^2\) can be calculated for each correlation coefficients to indicate the degree of heterogeneity.
  • Notes on the heterogeneity matrix (\(T^2\)) in stage 1 analysis:
    • The no. of effect sizes is usually much larger in MASEM than in conventional multivariate meta-analysis.
    • For example, it is uncommon to have \(p=5\) variables in a MASEM. There will be \(5 \times (5-1)/2=10\) effect sizes and \(10 \times (10+1)/2=55\) elements in \(T^2\).
    • There may not be enough studies to estimate the full \(T^2\).
    • If this happens, we may use a diagonal matrix for \(T^2\) (RE.type="Diag").

Random-effects TSSEM: Stage 2 analysis

  • Suppose the proposed correlation structure is \(\rho_\mathrm{R}(\theta)=vechs(P_\mathrm{R}(\theta))\), the discrepancy function is
    • \(F_\mathrm{WLS}(\theta) = (r_\mathrm{R} - \rho_\mathrm{R}(\theta))^T V_\mathrm{R}^{-1} (r_\mathrm{R} - \rho_\mathrm{R}(\theta))\)
    • The model is similar to that for the fixed-effects TSSEM.
  • Notes on the heterogeneity matrix (\(T^2\)) in stage 2 analysis:
    • \(T^2\) is not directly involved in the above fit function;
    • As \(V_\mathrm{R}\) is estimated after controlling for the random effects \(T^2\), \(V_\mathrm{R}\) has already been taken the random effects into account;
    • \(V_\mathrm{R}\) is usually larger than \(V_\mathrm{F}\). Thus, the SEs in the random-effects model are generally larger than those in the fixed-effects model.

OSMASEM approach

  • The current approaches, except the Bayesian approach, do not handle continuous moderators properly.
  • Jak and Cheung (2020) proposed a one-stage MASEM approach that handles categorical and continuous moderators in MASEM.
  • Suppose we want to use Lag to predict the path coefficient w2w, the model is
    • \(\beta_{w2w}\) = w2w_0 + w2w_1*Lag, where w2w_0 is the path coefficient when Lag=0 and w2w_1 is the change when Lag increases one unit.
Jak and Cheung (2020)
Jak and Cheung (2020)

Examples and applications

An example: A Five-factor model

  • This example illustrates how to apply MASEM to models with latent variables.
  • Digman (1997)30 reported 14 correlation matrices among the five-factor model.
    • Alpha factor (general personality of socialization): agreeableness (A), conscientiousness (C), and emotional stability (ES; the reverse of neuroticism)
    • Beta factor (personal growth): extraversion (E) and intellect (I)

The dataset

  • The dataset is available in the metaSEM package as Digman97 (see ?Digman97 after loading the metaSEM package).
  • To illustrate how to prepare data for MASEM, I have converted it to Excel format.
  • The first sheet (Info) includes the study characteristics, whereas the other sheets (Study 1 to Study 14) store the correlation matrices.

Reading the Excel file into R

## Load the library for MASEM
library(metaSEM)

## Load the library to read XLSX file
library(readxl)

## Read the study characteristics
study <- read_xlsx("Digman97.xlsx", sheet="Info")

## Display a few studies
head(study)
## # A tibble: 6 × 3
##   Study                               n Cluster     
##   <chr>                           <dbl> <chr>       
## 1 Digman 1 (1994)                   102 Children    
## 2 Digman 2 (1994)                   149 Children    
## 3 Digman 3 (1963c)                  334 Children    
## 4 Digman & Takemoto-Chock (1981b)   162 Children    
## 5 Graziano & Ward (1992)             91 Adolescents 
## 6 Yik & Bond (1993)                 656 Young adults
## Create an empty list to store the correlation matrices
Digman97.data <- list()
  
## Read 1 to 14 correlation matrices
for (i in 1:14) {
  ## Read each sheet and convert it into a matrix
  mat <- as.matrix(read_xlsx("Digman97.xlsx", sheet=paste0("Study ", i)))
  ## Add the row names
  rownames(mat) <- colnames(mat)
  ## Save it into a list
  Digman97.data[[i]] <- mat
}

## Add the names of the studies
names(Digman97.data) <- study$Study

## Show the first few studies
head(Digman97.data)
## $`Digman 1 (1994)`
##        A     C   ES     E    I
## A   1.00  0.62 0.41 -0.48 0.00
## C   0.62  1.00 0.59 -0.10 0.35
## ES  0.41  0.59 1.00  0.27 0.41
## E  -0.48 -0.10 0.27  1.00 0.37
## I   0.00  0.35 0.41  0.37 1.00
## 
## $`Digman 2 (1994)`
##        A    C   ES     E     I
## A   1.00 0.39 0.53 -0.30 -0.05
## C   0.39 1.00 0.59  0.07  0.44
## ES  0.53 0.59 1.00  0.09  0.22
## E  -0.30 0.07 0.09  1.00  0.45
## I  -0.05 0.44 0.22  0.45  1.00
## 
## $`Digman 3 (1963c)`
##       A     C   ES     E    I
## A  1.00  0.65 0.35  0.25 0.14
## C  0.65  1.00 0.37 -0.10 0.33
## ES 0.35  0.37 1.00  0.24 0.41
## E  0.25 -0.10 0.24  1.00 0.41
## I  0.14  0.33 0.41  0.41 1.00
## 
## $`Digman & Takemoto-Chock (1981b)`
##        A     C   ES     E     I
## A   1.00  0.65 0.70 -0.26 -0.03
## C   0.65  1.00 0.71 -0.16  0.24
## ES  0.70  0.71 1.00  0.01  0.11
## E  -0.26 -0.16 0.01  1.00  0.66
## I  -0.03  0.24 0.11  0.66  1.00
## 
## $`Graziano & Ward (1992)`
##       A    C   ES    E    I
## A  1.00 0.64 0.35 0.29 0.22
## C  0.64 1.00 0.27 0.16 0.22
## ES 0.35 0.27 1.00 0.32 0.36
## E  0.29 0.16 0.32 1.00 0.53
## I  0.22 0.22 0.36 0.53 1.00
## 
## $`Yik & Bond (1993)`
##       A    C   ES    E    I
## A  1.00 0.66 0.57 0.35 0.38
## C  0.66 1.00 0.45 0.20 0.31
## ES 0.57 0.45 1.00 0.49 0.31
## E  0.35 0.20 0.49 1.00 0.59
## I  0.38 0.31 0.31 0.59 1.00
## Extract the sample sizes
Digman97.n <- study$n
Digman97.n
##  [1]  102  149  334  162   91  656   70   70  277  227 1000  227   91 1040
## Extract the cluster
Digman97.cluster <- study$Cluster
Digman97.cluster
##  [1] "Children"      "Children"      "Children"      "Children"     
##  [5] "Adolescents"   "Young adults"  "Young adults"  "Young adults" 
##  [9] "Mature adults" "Mature adults" "Mature adults" "Mature adults"
## [13] "Mature adults" "Mature adults"
## Display the no. of studies
pattern.na(Digman97.data, show.na=FALSE)
##     A  C ES  E  I
## A  14 14 14 14 14
## C  14 14 14 14 14
## ES 14 14 14 14 14
## E  14 14 14 14 14
## I  14 14 14 14 14
## Display the cumulative sample sizes
pattern.n(Digman97.data, Digman97.n)
##       A    C   ES    E    I
## A  4496 4496 4496 4496 4496
## C  4496 4496 4496 4496 4496
## ES 4496 4496 4496 4496 4496
## E  4496 4496 4496 4496 4496
## I  4496 4496 4496 4496 4496

Proposed model

  • Most users find it easier to specify the models in a subset of the lavaan syntax.
  • The metaSEM package uses a RAM specification. We may convert the lavaan model to RAM specification with the lavaan2RAM() function.
model <- "## Factor loadings
          ## Alpha is measured by A, C, and ES
          Alpha =~ A + C + ES
          ## Beta is measured by E and I
          Beta =~ E + I
          ## Factor correlation between Alpha and Beta
          Alpha ~~ Beta"

## Display the model
plot(model, color="yellow")

## Convert the lavaan syntax into a RAM model as the metaSEM only knows the RAM model
## It is important to ensure that the variables are arranged in A, C, ES, E, and I.
RAM <- lavaan2RAM(model, obs.variables=c("A","C","ES","E","I"))
RAM
## $A
##       A   C   ES  E   I   Alpha           Beta         
## A     "0" "0" "0" "0" "0" "0.1*AONAlpha"  "0"          
## C     "0" "0" "0" "0" "0" "0.1*CONAlpha"  "0"          
## ES    "0" "0" "0" "0" "0" "0.1*ESONAlpha" "0"          
## E     "0" "0" "0" "0" "0" "0"             "0.1*EONBeta"
## I     "0" "0" "0" "0" "0" "0"             "0.1*IONBeta"
## Alpha "0" "0" "0" "0" "0" "0"             "0"          
## Beta  "0" "0" "0" "0" "0" "0"             "0"          
## 
## $S
##       A            C            ES             E            I           
## A     "0.5*AWITHA" "0"          "0"            "0"          "0"         
## C     "0"          "0.5*CWITHC" "0"            "0"          "0"         
## ES    "0"          "0"          "0.5*ESWITHES" "0"          "0"         
## E     "0"          "0"          "0"            "0.5*EWITHE" "0"         
## I     "0"          "0"          "0"            "0"          "0.5*IWITHI"
## Alpha "0"          "0"          "0"            "0"          "0"         
## Beta  "0"          "0"          "0"            "0"          "0"         
##       Alpha             Beta             
## A     "0"               "0"              
## C     "0"               "0"              
## ES    "0"               "0"              
## E     "0"               "0"              
## I     "0"               "0"              
## Alpha "1"               "0*AlphaWITHBeta"
## Beta  "0*AlphaWITHBeta" "1"              
## 
## $F
##    A C ES E I Alpha Beta
## A  1 0  0 0 0     0    0
## C  0 1  0 0 0     0    0
## ES 0 0  1 0 0     0    0
## E  0 0  0 1 0     0    0
## I  0 0  0 0 1     0    0
## 
## $M
##   A C ES E I Alpha Beta
## 1 0 0  0 0 0     0    0

Fixed-effects TSSEM

Stage 1 Analysis

  • The test statistic on testing the homogeneity test of the correlation matrix is \(\chi^2(df=130, N=4496)=1505.44, p < 0.001\).
  • Both the RMSEA and SRMR are very large, indicating that the proposed model (homogeneity of correlation matrices) did not fit the data well.
  • The estimates of the \(S[i,j]\) represent the element of the pooled correlation (or covariance) matrix.
## method="FEM": fixed-effects TSSEM
fixed1 <- tssem1(Digman97.data, Digman97.n, method="FEM")

## summary of the findings
summary(fixed1)
## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.363278  0.013368 27.1760 < 2.2e-16 ***
## S[1,3] 0.390375  0.012880 30.3077 < 2.2e-16 ***
## S[1,4] 0.103572  0.015047  6.8830 5.861e-12 ***
## S[1,5] 0.092286  0.015047  6.1331 8.621e-10 ***
## S[2,3] 0.416070  0.012519 33.2345 < 2.2e-16 ***
## S[2,4] 0.135148  0.014776  9.1464 < 2.2e-16 ***
## S[2,5] 0.141431  0.014866  9.5135 < 2.2e-16 ***
## S[3,4] 0.244459  0.014153 17.2724 < 2.2e-16 ***
## S[3,5] 0.138339  0.014834  9.3259 < 2.2e-16 ***
## S[4,5] 0.424566  0.012375 34.3071 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      4496.0000
## Chi-square of target model       1505.4443
## DF of target model                130.0000
## p value of target model             0.0000
## Chi-square of independence model 4471.4242
## DF of independence model          140.0000
## RMSEA                               0.1815
## RMSEA lower 95% CI                  0.1736
## RMSEA upper 95% CI                  0.1901
## SRMR                                0.1621
## TLI                                 0.6580
## CFI                                 0.6824
## AIC                              1245.4443
## BIC                               412.0217
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## extract coefficients
coef(fixed1)
##             A         C        ES         E          I
## A  1.00000000 0.3632782 0.3903748 0.1035716 0.09228557
## C  0.36327824 1.0000000 0.4160695 0.1351477 0.14143058
## ES 0.39037483 0.4160695 1.0000000 0.2444593 0.13833895
## E  0.10357155 0.1351477 0.2444593 1.0000000 0.42456626
## I  0.09228557 0.1414306 0.1383390 0.4245663 1.00000000

Stage 2 Analysis

  • Although the random-effects model is preferred, we illustrate how to fit the stage 2 analysis based on the fixed-effects model.
  • Important notes:
    • The standard errors (SEs) of the parameter estimates are underestimated (incorrect) when the correlation matrices are not homogeneous.
    • Readers should be cautious when interpreting the results of the stage 2 analysis when the homogeneity of the correlation matrices is rejected in the stage 1 analysis.
fixed2 <- tssem2(fixed1, RAM=RAM)
summary(fixed2)
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AONAlpha      0.562800  0.015380 0.532655 0.592944  36.593 < 2.2e-16 ***
## CONAlpha      0.605217  0.015324 0.575183 0.635251  39.495 < 2.2e-16 ***
## EONBeta       0.781413  0.034244 0.714296 0.848529  22.819 < 2.2e-16 ***
## ESONAlpha     0.719201  0.015685 0.688458 0.749943  45.852 < 2.2e-16 ***
## IONBeta       0.551374  0.026031 0.500354 0.602394  21.181 < 2.2e-16 ***
## AlphaWITHBeta 0.362678  0.022384 0.318806 0.406549  16.203 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                   65.4526
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           3112.7591
## DF of independence model                     10.0000
## RMSEA                                         0.0585
## RMSEA lower 95% CI                            0.0465
## RMSEA upper 95% CI                            0.0713
## SRMR                                          0.0284
## TLI                                           0.9505
## CFI                                           0.9802
## AIC                                          57.4526
## BIC                                          31.8088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
plot(fixed2, color="green")

Fixed-effects TSSEM with the type of participants as a moderator

Stage 1 analysis

  • As the assumption of the homogeneity of the correlation matrices has not been met, we group the studies into clusters based on the study characteristics.
  • For ease of illustration, we recode the data into younger and older participants.
  • From the results, splitting the data into younger and older participants does not explain the heterogeneity of the data well.
# Display the original study characteristic
table(Digman97.cluster)     
## Digman97.cluster
##   Adolescents      Children Mature adults  Young adults 
##             1             4             6             3
## Younger participants: "Children" and "Adolescents"
## Older participants: "Mature adults"
sample <- ifelse(Digman97.cluster %in% c("Children", "Adolescents"), 
                 yes="Younger participants", no="Older participants")

table(sample)
## sample
##   Older participants Younger participants 
##                    9                    5
## cluster: variable for the analysis with cluster
fixed1.cluster <- tssem1(Digman97.data, Digman97.n, method="FEM", cluster=sample)

summary(fixed1.cluster)
## $`Older participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##        Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.297471  0.015436 19.2716 < 2.2e-16 ***
## S[1,3] 0.370248  0.014532 25.4774 < 2.2e-16 ***
## S[1,4] 0.137694  0.016403  8.3944 < 2.2e-16 ***
## S[1,5] 0.098061  0.016724  5.8637 4.528e-09 ***
## S[2,3] 0.393692  0.014146 27.8306 < 2.2e-16 ***
## S[2,4] 0.183045  0.016055 11.4009 < 2.2e-16 ***
## S[2,5] 0.092774  0.016643  5.5743 2.485e-08 ***
## S[3,4] 0.260753  0.015554 16.7645 < 2.2e-16 ***
## S[3,5] 0.096141  0.016573  5.8009 6.597e-09 ***
## S[4,5] 0.411756  0.013900 29.6224 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                      3658.0000
## Chi-square of target model        825.9826
## DF of target model                 80.0000
## p value of target model             0.0000
## Chi-square of independence model 3000.9731
## DF of independence model           90.0000
## RMSEA                               0.1515
## RMSEA lower 95% CI                  0.1424
## RMSEA upper 95% CI                  0.1611
## SRMR                                0.1459
## TLI                                 0.7117
## CFI                                 0.7437
## AIC                               665.9826
## BIC                               169.6088
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## tssem1FEM(Cov = data.cluster[[i]], n = n.cluster[[i]], cor.analysis = cor.analysis, 
##     model.name = model.name, suppressWarnings = suppressWarnings)
## 
## Coefficients:
##         Estimate Std.Error z value  Pr(>|z|)    
## S[1,2]  0.604330  0.022125 27.3142 < 2.2e-16 ***
## S[1,3]  0.465536  0.027493 16.9327 < 2.2e-16 ***
## S[1,4] -0.031381  0.035940 -0.8732   0.38258    
## S[1,5]  0.061508  0.034547  1.7804   0.07500 .  
## S[2,3]  0.501432  0.026348 19.0311 < 2.2e-16 ***
## S[2,4] -0.060678  0.034557 -1.7559   0.07911 .  
## S[2,5]  0.320987  0.031064 10.3330 < 2.2e-16 ***
## S[3,4]  0.175437  0.033675  5.2097 1.891e-07 ***
## S[3,5]  0.305149  0.031586  9.6609 < 2.2e-16 ***
## S[4,5]  0.478640  0.026883 17.8045 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                      Value
## Sample size                       838.0000
## Chi-square of target model        346.2810
## DF of target model                 40.0000
## p value of target model             0.0000
## Chi-square of independence model 1470.4511
## DF of independence model           50.0000
## RMSEA                               0.2139
## RMSEA lower 95% CI                  0.1939
## RMSEA upper 95% CI                  0.2355
## SRMR                                0.1411
## TLI                                 0.7305
## CFI                                 0.7844
## AIC                               266.2810
## BIC                                77.0402
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)

Stage 2 analysis

  • As an illustration, we show how to fit the stage 2 analysis.
  • We should not trust the results because the estimates, e.g., Beta_I=3.28, for the younger participants are problematic.
fixed2.cluster <- tssem2(fixed1.cluster, RAM=RAM)

summary(fixed2.cluster)
## $`Older participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AONAlpha      0.512656  0.018206 0.476973 0.548340  28.158 < 2.2e-16 ***
## CONAlpha      0.549967  0.017945 0.514795 0.585140  30.647 < 2.2e-16 ***
## EONBeta       0.967136  0.058751 0.851986 1.082287  16.462 < 2.2e-16 ***
## ESONAlpha     0.732174  0.018929 0.695074 0.769273  38.680 < 2.2e-16 ***
## IONBeta       0.430649  0.029634 0.372568 0.488730  14.532 < 2.2e-16 ***
## AlphaWITHBeta 0.349236  0.028118 0.294125 0.404346  12.420 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                3658.0000
## Chi-square of target model                   21.9954
## DF of target model                            4.0000
## p value of target model                       0.0002
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2273.3179
## DF of independence model                     10.0000
## RMSEA                                         0.0351
## RMSEA lower 95% CI                            0.0217
## RMSEA upper 95% CI                            0.0500
## SRMR                                          0.0160
## TLI                                           0.9801
## CFI                                           0.9920
## AIC                                          13.9954
## BIC                                         -10.8233
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## 
## $`Younger participants`
## 
## Call:
## wls(Cov = coef.tssem1FEM(tssem1.obj), aCov = vcov.tssem1FEM(tssem1.obj), 
##     n = sum(tssem1.obj$n), RAM = RAM, Amatrix = Amatrix, Smatrix = Smatrix, 
##     Fmatrix = Fmatrix, diag.constraints = diag.constraints, cor.analysis = tssem1.obj$cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##                Estimate Std.Error    lbound    ubound z value Pr(>|z|)    
## AONAlpha       0.747647  0.023880  0.700842  0.794451 31.3081   <2e-16 ***
## CONAlpha       0.911705  0.019864  0.872772  0.950638 45.8969   <2e-16 ***
## EONBeta        0.152563  0.159128 -0.159322  0.464448  0.9587   0.3377    
## ESONAlpha      0.677435  0.025864  0.626743  0.728126 26.1926   <2e-16 ***
## IONBeta        3.283839  3.363262 -3.308033  9.875711  0.9764   0.3289    
## AlphaWITHBeta  0.117257  0.125421 -0.128565  0.363078  0.9349   0.3498    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                 838.0000
## Chi-square of target model                  145.6167
## DF of target model                            4.0000
## p value of target model                       0.0000
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model           2480.2403
## DF of independence model                     10.0000
## RMSEA                                         0.2057
## RMSEA lower 95% CI                            0.1778
## RMSEA upper 95% CI                            0.2350
## SRMR                                          0.1051
## TLI                                           0.8567
## CFI                                           0.9427
## AIC                                         137.6167
## BIC                                         118.6926
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Setup two plots side-by-side
layout(t(1:2))

## Plot the first group
plot(fixed2.cluster[[1]], col="green")
title("Younger participants")

## Plot the second group
plot(fixed2.cluster[[2]], col="yellow")
title("Older participants")

Random-effects TSSEM

Stage 1 Analysis

  • There are ten correlation coefficients in the analysis with 55 elements in the variance component matrix of the random effects.

  • We may specify RE.type="Diag" with only 10 variances in the diagonal variance component matrix.

  • The minimum and the maximum \(I^2\) are 0.7669 and 0.9326, respectively.

  • The random-effects model is more appropriate.

## method="REM": Random-effects model
random1 <- tssem1(Digman97.data, Digman97.n, method="REM", RE.type="Diag")

summary(random1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation (robust=FALSE)
## Coefficients:
##                Estimate   Std.Error      lbound      ubound z value  Pr(>|z|)
## Intercept1   0.38971908  0.05429256  0.28330762  0.49613054  7.1781 7.068e-13
## Intercept2   0.43265881  0.04145136  0.35141563  0.51390198 10.4377 < 2.2e-16
## Intercept3   0.04945631  0.06071079 -0.06953466  0.16844728  0.8146   0.41529
## Intercept4   0.09603708  0.04478711  0.00825595  0.18381822  2.1443   0.03201
## Intercept5   0.42724244  0.03911620  0.35057609  0.50390878 10.9224 < 2.2e-16
## Intercept6   0.11929318  0.04106203  0.03881309  0.19977328  2.9052   0.00367
## Intercept7   0.19292424  0.04757962  0.09966990  0.28617858  4.0548 5.018e-05
## Intercept8   0.22690164  0.03240892  0.16338132  0.29042196  7.0012 2.538e-12
## Intercept9   0.18105567  0.04258855  0.09758363  0.26452770  4.2513 2.126e-05
## Intercept10  0.43614968  0.03205960  0.37331402  0.49898535 13.6043 < 2.2e-16
## Tau2_1_1     0.03648372  0.01513055  0.00682839  0.06613905  2.4113   0.01590
## Tau2_2_2     0.01963098  0.00859789  0.00277942  0.03648253  2.2832   0.02242
## Tau2_3_3     0.04571438  0.01952285  0.00745030  0.08397846  2.3416   0.01920
## Tau2_4_4     0.02236122  0.00995083  0.00285794  0.04186449  2.2472   0.02463
## Tau2_5_5     0.01729072  0.00796404  0.00168149  0.03289995  2.1711   0.02992
## Tau2_6_6     0.01815481  0.00895896  0.00059557  0.03571405  2.0264   0.04272
## Tau2_7_7     0.02604881  0.01130266  0.00389602  0.04820161  2.3047   0.02119
## Tau2_8_8     0.00988761  0.00513713 -0.00018098  0.01995619  1.9247   0.05426
## Tau2_9_9     0.01988244  0.00895053  0.00233973  0.03742515  2.2214   0.02633
## Tau2_10_10   0.01049222  0.00501578  0.00066148  0.02032296  2.0918   0.03645
##                
## Intercept1  ***
## Intercept2  ***
## Intercept3     
## Intercept4  *  
## Intercept5  ***
## Intercept6  ** 
## Intercept7  ***
## Intercept8  ***
## Intercept9  ***
## Intercept10 ***
## Tau2_1_1    *  
## Tau2_2_2    *  
## Tau2_3_3    *  
## Tau2_4_4    *  
## Tau2_5_5    *  
## Tau2_6_6    *  
## Tau2_7_7    *  
## Tau2_8_8    .  
## Tau2_9_9    *  
## Tau2_10_10  *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 1220.334
## Degrees of freedom of the Q statistic: 130
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                               Estimate
## Intercept1: I2 (Q statistic)    0.9326
## Intercept2: I2 (Q statistic)    0.8872
## Intercept3: I2 (Q statistic)    0.9325
## Intercept4: I2 (Q statistic)    0.8703
## Intercept5: I2 (Q statistic)    0.8797
## Intercept6: I2 (Q statistic)    0.8480
## Intercept7: I2 (Q statistic)    0.8887
## Intercept8: I2 (Q statistic)    0.7669
## Intercept9: I2 (Q statistic)    0.8590
## Intercept10: I2 (Q statistic)   0.8193
## 
## Number of studies (or clusters): 14
## Number of observed statistics: 140
## Number of estimated parameters: 20
## Degrees of freedom: 120
## -2 log likelihood: -112.4196 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))
##  Intercept1  Intercept2  Intercept3  Intercept4  Intercept5  Intercept6 
##  0.38971908  0.43265881  0.04945631  0.09603708  0.42724244  0.11929318 
##  Intercept7  Intercept8  Intercept9 Intercept10 
##  0.19292424  0.22690164  0.18105567  0.43614968
## Average correlation matrix
## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
averageR <- vec2symMat(est_fixed, diag=FALSE)
dimnames(averageR) <- list(c("A", "C", "ES", "E", "I"),
                           c("A", "C", "ES", "E", "I"))
averageR
##             A         C        ES          E          I
## A  1.00000000 0.3897191 0.4326588 0.04945631 0.09603708
## C  0.38971908 1.0000000 0.4272424 0.11929318 0.19292424
## ES 0.43265881 0.4272424 1.0000000 0.22690164 0.18105567
## E  0.04945631 0.1192932 0.2269016 1.00000000 0.43614968
## I  0.09603708 0.1929242 0.1810557 0.43614968 1.00000000

Stage 2 Analysis

  • The stage 2 analysis is similar to that for the fixed-effects model.

  • We may trust the results based on the random-effects model.

  • The proposed model fits the data well with \(\chi^2(df=4, N=4496)=7.82, p=0.0984\), CFI=\(0.9922\), RMSEA=\(0.0146\), and SRMR=\(0.0436\).

  • The estimated correlation between the two factors (labeled AlphawithBeta in the printout) is \(0.3777\).

random2 <- tssem2(random1, RAM=RAM)
summary(random2)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##               Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## AONAlpha      0.569435  0.052425 0.466684 0.672187 10.8619 < 2.2e-16 ***
## CONAlpha      0.590630  0.052649 0.487439 0.693821 11.2182 < 2.2e-16 ***
## EONBeta       0.679955  0.075723 0.531541 0.828370  8.9795 < 2.2e-16 ***
## ESONAlpha     0.760455  0.061963 0.639009 0.881901 12.2726 < 2.2e-16 ***
## IONBeta       0.641842  0.072459 0.499825 0.783859  8.8580 < 2.2e-16 ***
## AlphaWITHBeta 0.377691  0.047402 0.284785 0.470596  7.9679 1.554e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                Value
## Sample size                                4496.0000
## Chi-square of target model                    7.8204
## DF of target model                            4.0000
## p value of target model                       0.0984
## Number of constraints imposed on "Smatrix"    0.0000
## DF manually adjusted                          0.0000
## Chi-square of independence model            501.6769
## DF of independence model                     10.0000
## RMSEA                                         0.0146
## RMSEA lower 95% CI                            0.0000
## RMSEA upper 95% CI                            0.0297
## SRMR                                          0.0436
## TLI                                           0.9806
## CFI                                           0.9922
## AIC                                          -0.1796
## BIC                                         -25.8234
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Plot the parameter estimates
plot(random2, color="green")

An example: A path model

  • Nohe et al. (2015)31 meta-analyzed the relationship between work-family conflict and strain.
  • The following shows a cross-lagged panel model on work (W) and strain (S).

Data preparation

  • The dataset is available as Nohe15A1 in the metaSEM package (see ?Nohe15A1 for the documentation). To illustrate how to read it as an external file, we have exported it as an Excel file Nohe15.xlsx.
library(readxl)
library(metaSEM)

## Read the Excel file
df <- read_xlsx("Nohe15.xlsx", sheet="Table A1")

head(df)
## # A tibble: 6 × 18
##   Study          N Participants Country Strain RelW1 RelW2 RelS1 RelS2 `% Women`
##   <chr>      <dbl> <chr>        <chr>   <chr>  <dbl> <dbl> <dbl> <dbl>     <dbl>
## 1 Britt & D…   489 Soldiers     United… Depre…  0.94  0.94  0.8   0.8         15
## 2 Demerouti…   335 Employment … the Ne… Exhau…  0.79  0.81  0.85  0.89        70
## 3 Ford (201…   328 Heterogeneo… United… Depre…  0.8   0.81  0.82  0.84        49
## 4 Hammer et…   234 Wives from … United… Depre…  0.91  0.91  0.9   0.9        100
## 5 Hammer et…   234 Husbands fr… United… Depre…  0.9   0.9   0.87  0.87         0
## 6 Innstrand…  2235 Professiona… Norway  Exhau…  0.7   0.7   0.86  0.88        46
## # ℹ 8 more variables: `Publication Status` <chr>, Lag <dbl>, `W1-S2` <dbl>,
## #   `S1-W2` <dbl>, `W1-S1` <dbl>, `W2-S2` <dbl>, `W1-W2` <dbl>, `S1-S2` <dbl>
## Variable names
my.var <- c("W1", "S1", "W2", "S2")

## Split each row as a list
my.list <- split(df, 1:nrow(df))

## Select the correlation coefficients and convert it into a correlation matrix
my.cor <- lapply(my.list, 
                 function(x) {## Convert the correlations into a correlation matrix
                              mat <- vec2symMat(unlist(x[c("W1-S1","W1-W2","W1-S2","S1-W2","S1-S2","W2-S2")]), diag = FALSE)
                              ## Add the dimensions for ease of reference
                              dimnames(mat) <- list(my.var, my.var)
                              ## Return the correlation matrix
                              mat})

## Add the study names for ease of reference
names(my.cor) <- df$Study

head(my.cor)
## $`Britt & Dawson (2005)`
##      W1   S1   W2   S2
## W1 1.00 0.29 0.58 0.22
## S1 0.29 1.00 0.24 0.57
## W2 0.58 0.24 1.00 0.27
## S2 0.22 0.57 0.27 1.00
## 
## $`Demerouti et al. (2004)`
##      W1   S1   W2   S2
## W1 1.00 0.53 0.57 0.41
## S1 0.53 1.00 0.41 0.68
## W2 0.57 0.41 1.00 0.54
## S2 0.41 0.68 0.54 1.00
## 
## $`Ford (2010)`
##      W1   S1   W2   S2
## W1 1.00 0.35 0.75 0.32
## S1 0.35 1.00 0.26 0.74
## W2 0.75 0.26 1.00 0.30
## S2 0.32 0.74 0.30 1.00
## 
## $`Hammer et al. (2005), female subsample`
##      W1   S1   W2   S2
## W1 1.00 0.32 0.57 0.22
## S1 0.32 1.00 0.30 0.43
## W2 0.57 0.30 1.00 0.30
## S2 0.22 0.43 0.30 1.00
## 
## $`Hammer et al. (2005), male subsample`
##      W1   S1   W2   S2
## W1 1.00 0.19 0.54 0.17
## S1 0.19 1.00 0.21 0.60
## W2 0.54 0.21 1.00 0.30
## S2 0.17 0.60 0.30 1.00
## 
## $`Innstrand et al. (2008)`
##      W1   S1   W2   S2
## W1 1.00 0.42 0.63 0.31
## S1 0.42 1.00 0.30 0.62
## W2 0.63 0.30 1.00 0.44
## S2 0.31 0.62 0.44 1.00
## Sample sizes
df$N
##  [1]  489  335  328  234  234 2235   76   94  236  239  239  138  160  151  409
## [16]   78  256  260  600  462  215 1292  470  665  403  153  201  382  130  946
## [31]  730   66
## Display the no. of studies
pattern.na(my.cor, show.na=FALSE)
##    W1 S1 W2 S2
## W1 32 32 32 32
## S1 32 32 32 32
## W2 32 32 32 32
## S2 32 32 32 32
## Display the cumulative sample sizes
pattern.n(my.cor, df$N)
##       W1    S1    W2    S2
## W1 12906 12906 12906 12906
## S1 12906 12906 12906 12906
## W2 12906 12906 12906 12906
## S2 12906 12906 12906 12906

TSSEM

Stage 1 analysis

  • We use both a fixed-effects and random-effects model to meta-analyze the correlation matrices. Since the fixed-effects model does not fit the data, we rely on the random-effects model.
  • The interpretations are similar to those in a multivariate meta-analysis.
    • The likelihood ratio (LR) statistic in testing the homogeneity of correlation matrices is \(\chi^2(df=186)=1,524, p < .001\), suggesting that the fixed-effects model does not fit the data well.
    • When we apply a random-effects model, the \(I^2\) varies from .77 to .88.
## Fit the stage one model (fixed-effects model)
fix1 <- tssem1(my.cor, df$N, method="FEM")

summary(fix1)
## 
## Call:
## tssem1FEM(Cov = Cov, n = n, cor.analysis = cor.analysis, model.name = model.name, 
##     cluster = cluster, suppressWarnings = suppressWarnings, silent = silent, 
##     run = run)
## 
## Coefficients:
##         Estimate Std.Error z value  Pr(>|z|)    
## S[1,2] 0.4239505 0.0073295  57.841 < 2.2e-16 ***
## S[1,3] 0.6242536 0.0054290 114.984 < 2.2e-16 ***
## S[1,4] 0.3349734 0.0078888  42.462 < 2.2e-16 ***
## S[2,3] 0.3288615 0.0079183  41.532 < 2.2e-16 ***
## S[2,4] 0.6315353 0.0053669 117.673 < 2.2e-16 ***
## S[3,4] 0.4342405 0.0072612  59.803 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                       Value
## Sample size                      12906.0000
## Chi-square of target model        1524.9269
## DF of target model                 186.0000
## p value of target model              0.0000
## Chi-square of independence model 18047.6811
## DF of independence model           192.0000
## RMSEA                                0.1336
## RMSEA lower 95% CI                   0.1276
## RMSEA upper 95% CI                   0.1400
## SRMR                                 0.1096
## TLI                                  0.9226
## CFI                                  0.9250
## AIC                               1152.9269
## BIC                               -235.6464
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Fit the stage one model (random-effects model)
rand1 <- tssem1(my.cor, df$N, method="REM")

summary(rand1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation (robust=FALSE)
## Coefficients:
##             Estimate Std.Error    lbound    ubound z value  Pr(>|z|)    
## Intercept1 0.3804522 0.0225616 0.3362323 0.4246720 16.8628 < 2.2e-16 ***
## Intercept2 0.6051298 0.0180362 0.5697794 0.6404802 33.5508 < 2.2e-16 ***
## Intercept3 0.3032290 0.0178803 0.2681842 0.3382738 16.9588 < 2.2e-16 ***
## Intercept4 0.3036392 0.0178408 0.2686718 0.3386065 17.0194 < 2.2e-16 ***
## Intercept5 0.6166503 0.0166427 0.5840312 0.6492694 37.0523 < 2.2e-16 ***
## Intercept6 0.3954085 0.0216645 0.3529469 0.4378701 18.2515 < 2.2e-16 ***
## Tau2_1_1   0.0134777 0.0038704 0.0058919 0.0210635  3.4823 0.0004972 ***
## Tau2_2_2   0.0087592 0.0025260 0.0038083 0.0137102  3.4676 0.0005252 ***
## Tau2_3_3   0.0071123 0.0022470 0.0027082 0.0115163  3.1652 0.0015496 ** 
## Tau2_4_4   0.0070585 0.0022121 0.0027229 0.0113941  3.1909 0.0014183 ** 
## Tau2_5_5   0.0072634 0.0021092 0.0031293 0.0113974  3.4436 0.0005740 ***
## Tau2_6_6   0.0122813 0.0034848 0.0054513 0.0191114  3.5243 0.0004246 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 1466.163
## Degrees of freedom of the Q statistic: 186
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.8829
## Intercept2: I2 (Q statistic)   0.8973
## Intercept3: I2 (Q statistic)   0.7743
## Intercept4: I2 (Q statistic)   0.7718
## Intercept5: I2 (Q statistic)   0.8810
## Intercept6: I2 (Q statistic)   0.8748
## 
## Number of studies (or clusters): 32
## Number of observed statistics: 192
## Number of estimated parameters: 12
## Degrees of freedom: 180
## -2 log likelihood: -300.1701 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Extract the fixed-effects estimates
R1 <- coef(rand1, select = "fixed")
R1
## Intercept1 Intercept2 Intercept3 Intercept4 Intercept5 Intercept6 
##  0.3804522  0.6051298  0.3032290  0.3036392  0.6166503  0.3954085
## Convert it to a correlation matrix
R1 <- vec2symMat(R1, diag = FALSE)
R1
##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.3804522 0.6051298 0.3032290
## [2,] 0.3804522 1.0000000 0.3036392 0.6166503
## [3,] 0.6051298 0.3036392 1.0000000 0.3954085
## [4,] 0.3032290 0.6166503 0.3954085 1.0000000
## Add the dimension names for ease of reference
dimnames(R1) <- list(my.var, my.var)
R1
##           W1        S1        W2        S2
## W1 1.0000000 0.3804522 0.6051298 0.3032290
## S1 0.3804522 1.0000000 0.3036392 0.6166503
## W2 0.6051298 0.3036392 1.0000000 0.3954085
## S2 0.3032290 0.6166503 0.3954085 1.0000000

Stage 2 analysis

  • We then test the cross-lagged model.
  • As the proposed model fits the data perfectly, we ignore the fit indices in the output.
## Model in lavaan syntax
## The parameter labels are not necessary. But they make the output easier to follow.
model1 <- "W2 ~ w2w*W1 + s2w*S1    # path coefficeients   
           S2 ~ s2s*S1 + w2s*W1
           S1 ~~ w1withs1*W1       # correlation
           S2 ~~ w2withs2*W2
           W1 ~~ 1*W1              # variances of the independent variables are fixed at 1
           S1 ~~ 1*S1"

## Display the model
## See the layout argument in ?semPaths
## layout can be tree, circle, spring, tree2, circle2
plot(model1, color="yellow", layout="spring")

## Convert the above model to RAM specification used by metaSEM.
## We also have to specify how the variables are arranged in the data.
RAM1 <- lavaan2RAM(model1, obs.variables = my.var)
RAM1
## $A
##    W1        S1        W2  S2 
## W1 "0"       "0"       "0" "0"
## S1 "0"       "0"       "0" "0"
## W2 "0.1*w2w" "0.1*s2w" "0" "0"
## S2 "0.1*w2s" "0.1*s2s" "0" "0"
## 
## $S
##    W1           S1           W2             S2            
## W1 "1"          "0*w1withs1" "0"            "0"           
## S1 "0*w1withs1" "1"          "0"            "0"           
## W2 "0"          "0"          "0.5*W2WITHW2" "0*w2withs2"  
## S2 "0"          "0"          "0*w2withs2"   "0.5*S2WITHS2"
## 
## $F
##    W1 S1 W2 S2
## W1  1  0  0  0
## S1  0  1  0  0
## W2  0  0  1  0
## S2  0  0  0  1
## 
## $M
##   W1 S1 W2 S2
## 1  0  0  0  0
## Fit the stage two model
rand2a <- tssem2(rand1, RAM=RAM1)

summary(rand2a)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##          Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## s2s      0.586124  0.020790 0.545376 0.626872 28.1926 < 2.2e-16 ***
## w2s      0.080237  0.024842 0.031547 0.128927  3.2299 0.0012385 ** 
## s2w      0.085841  0.024796 0.037242 0.134440  3.4619 0.0005364 ***
## w2w      0.572471  0.022265 0.528834 0.616109 25.7122 < 2.2e-16 ***
## w1withs1 0.380452  0.022562 0.336232 0.424672 16.8628 < 2.2e-16 ***
## w2withs2 0.168885  0.025232 0.119431 0.218338  6.6933 2.182e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                            Value
## Sample size                                12906
## Chi-square of target model                     0
## DF of target model                             0
## p value of target model                        0
## Number of constraints imposed on "Smatrix"     0
## DF manually adjusted                           0
## Chi-square of independence model            3079
## DF of independence model                       6
## RMSEA                                          0
## RMSEA lower 95% CI                             0
## RMSEA upper 95% CI                             0
## SRMR                                           0
## TLI                                         -Inf
## CFI                                            1
## AIC                                            0
## BIC                                            0
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
plot(rand2a, color="green", layout="spring")

Stage 2 analysis with constraints

  • Suppose we have reasons to test whether the same variable pair (s2s and w2w) are the same and also the different variable pair (s2w and w2s) are the same, we may test them by imposing equality constraints.
  • As these two models are nested, we may use an LR test to compare them. The LR statistic is \(\chi^2(df=2)=0.22, p=.89\), suggesting that there is not enough evidence to reject the null hypothesis that the same pair and different pair effects.
model2 <- "W2 ~ same*W1 + diff*S1     # regression coefficeients   
           S2 ~ same*S1 + diff*W1
           S1 ~~ w1cs1*W1             # correlation
           S2 ~~ w2cs2*W2
           W1 ~~ 1*W1                 # variances of the independent variables are fixed at 1
           S1 ~~ 1*S1"

## Display the model
plot(model2, color="yellow")

## Convert the above model to RAM specification used by metaSEM.
## We also have to specify how the variables are arranged in the data.
RAM2 <- lavaan2RAM(model2, obs.variables = my.var)

## Fit the stage two model
rand2b <- tssem2(rand1, RAM=RAM2)

summary(rand2b)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation
## Coefficients:
##       Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
## diff  0.082840  0.019712 0.044206 0.121474  4.2026 2.638e-05 ***
## same  0.579844  0.015217 0.550020 0.609669 38.1053 < 2.2e-16 ***
## w1cs1 0.380451  0.022562 0.336231 0.424671 16.8628 < 2.2e-16 ***
## w2cs2 0.168824  0.025241 0.119353 0.218295  6.6886 2.254e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Goodness-of-fit indices:
##                                                 Value
## Sample size                                12906.0000
## Chi-square of target model                     0.2247
## DF of target model                             2.0000
## p value of target model                        0.8937
## Number of constraints imposed on "Smatrix"     0.0000
## DF manually adjusted                           0.0000
## Chi-square of independence model            3078.9540
## DF of independence model                       6.0000
## RMSEA                                          0.0000
## RMSEA lower 95% CI                             0.0000
## RMSEA upper 95% CI                             0.0079
## SRMR                                           0.0033
## TLI                                            1.0017
## CFI                                            1.0000
## AIC                                           -3.7753
## BIC                                          -18.7062
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Comparing nested models with anova()
anova(rand2a, rand2b)
##                 base         comparison ep     minus2LL df       AIC    diffLL
## 1 TSSEM2 Correlation               <NA>  6 1.245611e-19 -6 12.000000        NA
## 2 TSSEM2 Correlation TSSEM2 Correlation  4 2.246694e-01 -4  8.224669 0.2246694
##   diffdf         p
## 1     NA        NA
## 2      2 0.8937451
plot(rand2b, color="green")

OSMASEM: Moderator analysis

  • We may also want to check if the covariates help us to explain the heterogeneity in the data.
  • Lag is the time lag between the measurement waves in months. It is reasonable to expect that it may be related to the strength of the paths.
    • By comparing the models with and without the covariate, the LR statistic is \(\chi^2(df=4)=23.52, p<.001\), suggesting that at least one coefficient is non-zero.
    • In general, the effect is negative (rows 7-10) in summary(fit1), indicating that the paths become weaker over time.
    • By inspecting rows 7-10 in summary(fit1), only w2w_1 (row 7) is statistically significant.
    • We may compute the simple slopes at the mean, mean $$1 SD of Lag. Since Lag is standardized, it makes it easier to calculate the effects. The path coefficients at -1v SD, mean, +1 SD are .64, .57, and .51, respectively.
## Convert the correlation matrices into a dataframe, which is required in osmasem()
Nohe.df <- Cor2DataFrame(x=my.cor, n=df$N)

## Standardize the moderator "Lag" to improve numerical stability and add it into the dataframe
Nohe.df$data$Lag <- scale(df$Lag)
    
head(Nohe.df$data)
##                                        S1_W1 W2_W1 S2_W1 W2_S1 S2_S1 S2_W2
## Britt...Dawson..2005.                   0.29  0.58  0.22  0.24  0.57  0.27
## Demerouti.et.al...2004.                 0.53  0.57  0.41  0.41  0.68  0.54
## Ford..2010.                             0.35  0.75  0.32  0.26  0.74  0.30
## Hammer.et.al...2005...female.subsample  0.32  0.57  0.22  0.30  0.43  0.30
## Hammer.et.al...2005...male.subsample    0.19  0.54  0.17  0.21  0.60  0.30
## Innstrand.et.al...2008.                 0.42  0.63  0.31  0.30  0.62  0.44
##                                        C(S1_W1 S1_W1) C(W2_W1 S1_W1)
## Britt...Dawson..2005.                     0.001422479   0.0002033635
## Demerouti.et.al...2004.                   0.002076395   0.0002968500
## Ford..2010.                               0.002120708   0.0003031852
## Hammer.et.al...2005...female.subsample    0.002972616   0.0004249776
## Hammer.et.al...2005...male.subsample      0.002972616   0.0004249776
## Innstrand.et.al...2008.                   0.000311227   0.0000444943
##                                        C(S2_W1 S1_W1) C(W2_S1 S1_W1)
## Britt...Dawson..2005.                    0.0008791243   0.0008736611
## Demerouti.et.al...2004.                  0.0012832591   0.0012752845
## Ford..2010.                              0.0013106457   0.0013025009
## Hammer.et.al...2005...female.subsample   0.0018371444   0.0018257278
## Hammer.et.al...2005...male.subsample     0.0018371444   0.0018257278
## Innstrand.et.al...2008.                  0.0001923453   0.0001911500
##                                        C(S2_S1 S1_W1) C(S2_W2 S1_W1)
## Britt...Dawson..2005.                    2.047346e-04   0.0004890894
## Demerouti.et.al...2004.                  2.988513e-04   0.0007139245
## Ford..2010.                              3.052293e-04   0.0007291607
## Hammer.et.al...2005...female.subsample   4.278427e-04   0.0010220714
## Hammer.et.al...2005...male.subsample     4.278427e-04   0.0010220714
## Innstrand.et.al...2008.                  4.479427e-05   0.0001070088
##                                        C(W2_W1 W2_W1) C(S2_W1 W2_W1)
## Britt...Dawson..2005.                    0.0007981106   3.750391e-04
## Demerouti.et.al...2004.                  0.0011650033   5.474451e-04
## Ford..2010.                              0.0011898662   5.591284e-04
## Hammer.et.al...2005...female.subsample   0.0016678466   7.837355e-04
## Hammer.et.al...2005...male.subsample     0.0016678466   7.837355e-04
## Innstrand.et.al...2008.                  0.0001746202   8.205553e-05
##                                        C(W2_S1 W2_W1) C(S2_S1 W2_W1)
## Britt...Dawson..2005.                    3.671018e-04   1.042810e-04
## Demerouti.et.al...2004.                  5.358591e-04   1.522191e-04
## Ford..2010.                              5.472951e-04   1.554677e-04
## Hammer.et.al...2005...female.subsample   7.671487e-04   2.179205e-04
## Hammer.et.al...2005...male.subsample     7.671487e-04   2.179205e-04
## Innstrand.et.al...2008.                  8.031892e-05   2.281583e-05
##                                        C(S2_W2 W2_W1) C(S2_W1 S2_W1)
## Britt...Dawson..2005.                    2.035379e-04    0.001649681
## Demerouti.et.al...2004.                  2.971046e-04    0.002408042
## Ford..2010.                              3.034452e-04    0.002459434
## Hammer.et.al...2005...female.subsample   4.253420e-04    0.003447411
## Hammer.et.al...2005...male.subsample     4.253420e-04    0.003447411
## Innstrand.et.al...2008.                  4.453246e-05    0.000360937
##                                        C(W2_S1 S2_W1) C(S2_S1 S2_W1)
## Britt...Dawson..2005.                    0.0005769098   3.592773e-04
## Demerouti.et.al...2004.                  0.0008421160   5.244375e-04
## Ford..2010.                              0.0008600880   5.356298e-04
## Hammer.et.al...2005...female.subsample   0.0012055935   7.507973e-04
## Hammer.et.al...2005...male.subsample     0.0012055935   7.507973e-04
## Innstrand.et.al...2008.                  0.0001262232   7.860697e-05
##                                        C(S2_W2 S2_W1) C(W2_S1 W2_S1)
## Britt...Dawson..2005.                    0.0008607786   0.0016601406
## Demerouti.et.al...2004.                  0.0012564798   0.0024233097
## Ford..2010.                              0.0012832949   0.0024750266
## Hammer.et.al...2005...female.subsample   0.0017988065   0.0034692681
## Hammer.et.al...2005...male.subsample     0.0017988065   0.0034692681
## Innstrand.et.al...2008.                  0.0001883314   0.0003632254
##                                        C(S2_S1 W2_S1) C(S2_W2 W2_S1)
## Britt...Dawson..2005.                    3.727412e-04   0.0008739022
## Demerouti.et.al...2004.                  5.440909e-04   0.0012756364
## Ford..2010.                              5.557026e-04   0.0013028604
## Hammer.et.al...2005...female.subsample   7.789335e-04   0.0018262316
## Hammer.et.al...2005...male.subsample     7.789335e-04   0.0018262316
## Innstrand.et.al...2008.                  8.155277e-05   0.0001912028
##                                        C(S2_S1 S2_S1) C(S2_W2 S2_S1)
## Britt...Dawson..2005.                    0.0007805638   0.0001951863
## Demerouti.et.al...2004.                  0.0011393902   0.0002849138
## Ford..2010.                              0.0011637064   0.0002909943
## Hammer.et.al...2005...female.subsample   0.0016311782   0.0004078894
## Hammer.et.al...2005...male.subsample     0.0016311782   0.0004078894
## Innstrand.et.al...2008.                  0.0001707811   0.0000427052
##                                        C(S2_W2 S2_W2)        Lag
## Britt...Dawson..2005.                    0.0013981148 -0.6794521
## Demerouti.et.al...2004.                  0.0020408303 -0.7711151
## Ford..2010.                              0.0020843846 -0.8016694
## Hammer.et.al...2005...female.subsample   0.0029217015 -0.1294740
## Hammer.et.al...2005...male.subsample     0.0029217015 -0.1294740
## Innstrand.et.al...2008.                  0.0003058963  0.6038301
## Fit the model without any moderator
## The results are similat to that in TSSEM.
fit0 <- osmasem(model.name="No moderator", RAM=RAM1, data=Nohe.df)
summary(fit0)
## Summary of No moderator 
##  
## free parameters:
##        name  matrix row col    Estimate  Std.Error A    z value     Pr(>|z|)
## 1       w2w      A0  W2  W1  0.57247133 0.02226456    25.712219 0.000000e+00
## 2       w2s      A0  S2  W1  0.08023683 0.02484214     3.229868 1.238476e-03
## 3       s2w      A0  W2  S1  0.08584122 0.02479590     3.461912 5.363533e-04
## 4       s2s      A0  S2  S1  0.58612403 0.02079000    28.192589 0.000000e+00
## 5  w1withs1      S0  S1  W1  0.38045217 0.02256156    16.862847 0.000000e+00
## 6  w2withs2      S0  S2  W2  0.16888459 0.02523192     6.693291 2.182055e-11
## 7    Tau1_1 vecTau1   1   1 -4.30671722 0.28716838   -14.997184 0.000000e+00
## 8    Tau1_2 vecTau1   2   1 -4.73764524 0.28838625   -16.428125 0.000000e+00
## 9    Tau1_3 vecTau1   3   1 -4.94593056 0.31593254   -15.655021 0.000000e+00
## 10   Tau1_4 vecTau1   4   1 -4.95352326 0.31339094   -15.806211 0.000000e+00
## 11   Tau1_5 vecTau1   5   1 -4.92491364 0.29039425   -16.959405 0.000000e+00
## 12   Tau1_6 vecTau1   6   1 -4.39967581 0.28374629   -15.505668 0.000000e+00
## 
## To obtain confidence intervals re-run with intervals=TRUE
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             12                    180             -300.1701
##    Saturated:             27                    165                    NA
## Independence:             12                    180                    NA
## Number of observations/statistics: 12906/192
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -660.1701              -276.1701                -276.1459
## BIC:     -2003.9507              -186.5847                -224.7195
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-07-21 15:01:14 
## Wall clock time: 0.0996995 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
plot(fit0)

## Create a matrix to present "Lag" moderator
## We use Lag to moderate the following paths:
## W1 -> W2
## W1 -> S2
## S1 -> W2
## S1 -> S2
A1 <- create.modMatrix(RAM=RAM1, output="A", mod="Lag")
A1
##    W1           S1           W2  S2 
## W1 "0"          "0"          "0" "0"
## S1 "0"          "0"          "0" "0"
## W2 "0*data.Lag" "0*data.Lag" "0" "0"
## S2 "0*data.Lag" "0*data.Lag" "0" "0"
## Fit the model with Lag as a moderator
fit1 <- osmasem(model.name="Lag as a moderator", RAM=RAM1, Ax=A1, data=Nohe.df)
summary(fit1)
## Summary of Lag as a moderator 
##  
## free parameters:
##        name  matrix row col     Estimate  Std.Error A     z value     Pr(>|z|)
## 1       w2w      A0  W2  W1  0.573039774 0.01839766    31.1474347 0.000000e+00
## 2       w2s      A0  S2  W1  0.079844963 0.02419489     3.3000756 9.665879e-04
## 3       s2w      A0  W2  S1  0.085391026 0.02393613     3.5674540 3.604667e-04
## 4       s2s      A0  S2  S1  0.586234249 0.01962561    29.8708837 0.000000e+00
## 5  w1withs1      S0  S1  W1  0.381183701 0.02282277    16.7019024 0.000000e+00
## 6  w2withs2      S0  S2  W2  0.166974817 0.02500439     6.6778197 2.425238e-11
## 7     w2w_1      A1  W2  W1 -0.062015587 0.01850574    -3.3511549 8.047527e-04
## 8     w2s_1      A1  S2  W1 -0.025933718 0.02096724    -1.2368685 2.161359e-01
## 9     s2w_1      A1  W2  S1 -0.002382818 0.02055897    -0.1159016 9.077305e-01
## 10    s2s_1      A1  S2  S1 -0.027809761 0.01974172    -1.4086798 1.589299e-01
## 11   Tau1_1 vecTau1   1   1 -4.276381062 0.28720205   -14.8898001 0.000000e+00
## 12   Tau1_2 vecTau1   2   1 -5.261036068 0.32310515   -16.2827364 0.000000e+00
## 13   Tau1_3 vecTau1   3   1 -5.048387485 0.32015075   -15.7687824 0.000000e+00
## 14   Tau1_4 vecTau1   4   1 -5.039816725 0.31967923   -15.7652303 0.000000e+00
## 15   Tau1_5 vecTau1   5   1 -5.074951110 0.29424551   -17.2473355 0.000000e+00
## 16   Tau1_6 vecTau1   6   1 -4.397726436 0.28288196   -15.5461536 0.000000e+00
## 
## To obtain confidence intervals re-run with intervals=TRUE
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             16                    176             -323.6921
##    Saturated:             27                    165                    NA
## Independence:             12                    180                    NA
## Number of observations/statistics: 12906/192
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:      -675.6921              -291.6921                -291.6499
## BIC:     -1989.6109              -172.2449                -223.0913
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2024-07-21 15:01:14 
## Wall clock time: 0.2052836 secs 
## optimizer:  SLSQP 
## OpenMx version number: 2.21.11 
## Need help?  See help(mxSummary)
## Comparing the models with and without the covariate
anova(fit1, fit0)
##                 base   comparison ep  minus2LL  df       AIC   diffLL diffdf
## 1 Lag as a moderator         <NA> 16 -323.6921 176 -291.6921       NA     NA
## 2 Lag as a moderator No moderator 12 -300.1701 180 -276.1701 23.52199      4
##              p
## 1           NA
## 2 9.957486e-05
## Effect of w2w when Lag is at the mean value
mxEval(w2w, fit1$mx.fit)
##           [,1]
## [1,] 0.5730398
## Effect of w2w when Lag is at +1 SD as Lag is standardized
mxEval(w2w + w2w_1, fit1$mx.fit)
##           [,1]
## [1,] 0.5110242
## Effect of w2w when Lag is at -1 SD as Lag is standardized
mxEval(w2w - w2w_1, fit1$mx.fit)
##           [,1]
## [1,] 0.6350554

An example: A mediation model

  • Schutte, Keng, and Cheung (2021)32 examined the association between dispositional mindfulness and gratitude, dispositional mindfulness and emotional intelligence, and emotional intelligence and gratitude via a mediation model. It involved 41 samples.
  • Although we may individually meta-analyze the bivariate relationships between EI and mindfulness, EI and gratitude, and mindfulness and gratitude, the results are not ideal. Researchers want to address the mediating role of the EI.
  • Another feature of this study is that there is no study including all three variables. Therefore, this study can address research questions beyond the primary studies. However, one limitation is that there may be a large degree of heterogeneity. Future studies may test the validity of the mediation model by including all three variables.

Data preparation

library(metaSEM)

## Read the SPSS dataset
my.df <- foreign::read.spss("Schutte21.sav", use.value.labels = TRUE, to.data.frame=TRUE)

## A function to convert rows into a 3x3 correlation matrix
create.matrix <- function(x, type=c(1, 2, 3)) {
  mat <- matrix(NA, ncol=3, nrow=3)
  diag(mat) <- 1
  type <- as.character(type)
  ## Mindfulness, EI, Gratitude
  ## 1: Mindfulness and EI
  ## 2: Mindfulness and Gratitude
  ## 3: EI and Gratitude
  switch(type,
         "1" = mat[1, 2] <- mat[2, 1] <- unlist(x),
         "2" = mat[1, 3] <- mat[3, 1] <- unlist(x),
         "3" = mat[2, 3] <- mat[3, 2] <- unlist(x))
  mat
}

## Convert the data to correlation matrices using the create.matrix().
my.cor <- lapply(split(my.df, seq(nrow(my.df))),
                 function(x, y) create.matrix(x["Effect_size"], x["Type_of_Association"]))

## Variable names
varlist <- c("Mindfulness", "EI", "Gratitude")

## Add the variable names in the dimnames.
my.cor <- lapply(my.cor, function(x) {dimnames(x) <- list(varlist, varlist); x}  )

## Add the study names
names(my.cor) <- my.df$Study_name

## Correlation matrices in the analysis
head(my.cor)
## $`Aras, 2015                                        `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.39        NA
## EI                 0.39 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Bao et al., 2015                                  `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.34        NA
## EI                 0.34 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Bravo et al., 2015                                `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.13        NA
## EI                 0.13 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Brown and Ryan, 2003_Sample A                     `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.46        NA
## EI                 0.46 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Brown and Ryan, 2003_Sample B                     `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.42        NA
## EI                 0.42 1.00        NA
## Gratitude            NA   NA         1
## 
## $`Brown and Ryan, 2003_Sample C                     `
##             Mindfulness   EI Gratitude
## Mindfulness        1.00 0.37        NA
## EI                 0.37 1.00        NA
## Gratitude            NA   NA         1
## Sample sizes
my.n <- my.df$N
my.n
##  [1]  100  380   83  313  327  207 1157  190  256  108  801  365   39   60  250
## [16]  427  470  228  375  234  121  535  156  970  123   99 1392  125  124  200
## [31]  200   72  700  341  313  103  406  299  321  327  200
## Number of studies in each cell
pattern.na(my.cor, show.na = FALSE)
##             Mindfulness EI Gratitude
## Mindfulness          41 26         8
## EI                   26 41         7
## Gratitude             8  7        41
## Total sample sizes in each cell
pattern.n(my.cor, my.n)
##             Mindfulness    EI Gratitude
## Mindfulness       13497  6369      3130
## EI                 6369 13497      3998
## Gratitude          3130  3998     13497

First-stage of analysis

## Stage 1 analysis: find an average correlation matrix
stage1 <- tssem1(my.cor, my.n)
summary(stage1)
## 
## Call:
## meta(y = ES, v = acovR, RE.constraints = Diag(paste0(RE.startvalues, 
##     "*Tau2_", 1:no.es, "_", 1:no.es)), RE.lbound = RE.lbound, 
##     I2 = I2, model.name = model.name, suppressWarnings = TRUE, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: z statistic approximation (robust=FALSE)
## Coefficients:
##              Estimate  Std.Error     lbound     ubound z value  Pr(>|z|)    
## Intercept1  0.4007062  0.0342472  0.3335830  0.4678294 11.7004 < 2.2e-16 ***
## Intercept2  0.2188135  0.0290679  0.1618414  0.2757855  7.5277 5.174e-14 ***
## Intercept3  0.3123689  0.0416955  0.2306472  0.3940906  7.4917 6.795e-14 ***
## Tau2_1_1    0.0260509  0.0089658  0.0084783  0.0436236  2.9056  0.003666 ** 
## Tau2_2_2    0.0038767  0.0030542 -0.0021094  0.0098628  1.2693  0.204327    
## Tau2_3_3    0.0095986  0.0058884 -0.0019424  0.0211396  1.6301  0.103083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Q statistic on the homogeneity of effect sizes: 262.1055
## Degrees of freedom of the Q statistic: 38
## P value of the Q statistic: 0
## 
## Heterogeneity indices (based on the estimated Tau2):
##                              Estimate
## Intercept1: I2 (Q statistic)   0.9209
## Intercept2: I2 (Q statistic)   0.5784
## Intercept3: I2 (Q statistic)   0.7998
## 
## Number of studies (or clusters): 41
## Number of observed statistics: 41
## Number of estimated parameters: 6
## Degrees of freedom: 35
## -2 log likelihood: -44.57614 
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values may indicate problems.)
## Average correlation matrix
meanR <- vec2symMat(coef(stage1, select = "fixed"), diag = FALSE)
dimnames(meanR) <- list(varlist, varlist)
meanR
##             Mindfulness        EI Gratitude
## Mindfulness   1.0000000 0.4007062 0.2188135
## EI            0.4007062 1.0000000 0.3123689
## Gratitude     0.2188135 0.3123689 1.0000000
## Absolute heterogeneity variance: tau^2
tau2 <- vec2symMat(coef(stage1, select = "random"), diag = FALSE)
dimnames(tau2) <- list(varlist, varlist)
tau2
##             Mindfulness          EI   Gratitude
## Mindfulness  1.00000000 0.026050950 0.003876740
## EI           0.02605095 1.000000000 0.009598596
## Gratitude    0.00387674 0.009598596 1.000000000
## Relative heterogeneity index: I^2
I2 <- vec2symMat(summary(stage1)$I2.values[, "Estimate"], diag = FALSE)
dimnames(I2) <- list(varlist, varlist)
I2
##             Mindfulness       EI Gratitude
## Mindfulness   1.0000000 0.920929 0.5783717
## EI            0.9209290 1.000000 0.7998260
## Gratitude     0.5783717 0.799826 1.0000000

Second-stage of analysis

## Proposed model
model <- "## cp: c prime (c'), the common notation for the direct effect.
          Gratitude ~ cp*Mindfulness + b*EI
          EI ~ a*Mindfulness
          Mindfulness ~~ 1*Mindfulness
          ## Define direct, indirect, and total effects
          Direct := cp
          Indirect := a*b
          Total := a*b + cp"

plot(model, color="yellow")

RAM1 <- lavaan2RAM(model, obs.variables = varlist)
RAM1
## $A
##             Mindfulness EI      Gratitude
## Mindfulness "0"         "0"     "0"      
## EI          "0.1*a"     "0"     "0"      
## Gratitude   "0.1*cp"    "0.1*b" "0"      
## 
## $S
##             Mindfulness EI             Gratitude                   
## Mindfulness "1"         "0"            "0"                         
## EI          "0"         "0.5*EIWITHEI" "0"                         
## Gratitude   "0"         "0"            "0.5*GratitudeWITHGratitude"
## 
## $F
##             Mindfulness EI Gratitude
## Mindfulness           1  0         0
## EI                    0  1         0
## Gratitude             0  0         1
## 
## $M
##   Mindfulness EI Gratitude
## 1           0  0         0
## 
## $mxalgebras
## $mxalgebras$Direct
## mxAlgebra 'Direct' 
## $formula:  cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Indirect
## mxAlgebra 'Indirect' 
## $formula:  a * b 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Total
## mxAlgebra 'Total' 
## $formula:  a * b + cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## Stage 2 analysis: fit the path model
## Likelihood-based CI: intervals.type = "LB"
stage2 <- tssem2(stage1, RAM=RAM1, intervals.type = "LB")
summary(stage2)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
##    Estimate Std.Error   lbound   ubound z value Pr(>|z|)
## a  0.400706        NA 0.333466 0.467842      NA       NA
## b  0.267667        NA 0.166085 0.369322      NA       NA
## cp 0.111558        NA 0.028737 0.190505      NA       NA
## 
## mxAlgebras objects (and their 95% likelihood-based CIs):
##                   lbound  Estimate    ubound
## Direct[1,1]   0.02873693 0.1115576 0.1905046
## Indirect[1,1] 0.06580411 0.1072559 0.1559166
## Total[1,1]    0.16178591 0.2188135 0.2760576
## 
## Goodness-of-fit indices:
##                                               Value
## Sample size                                13497.00
## Chi-square of target model                     0.00
## DF of target model                             0.00
## p value of target model                        0.00
## Number of constraints imposed on "Smatrix"     0.00
## DF manually adjusted                           0.00
## Chi-square of independence model             249.69
## DF of independence model                       3.00
## RMSEA                                          0.00
## RMSEA lower 95% CI                             0.00
## RMSEA upper 95% CI                             0.00
## SRMR                                           0.00
## TLI                                            -Inf
## CFI                                            1.00
## AIC                                            0.00
## BIC                                            0.00
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
plot(stage2, color="green")

Testing the direct effect equals the indirect effect (advanced topic!!!)

  • The previous analysis seems to indicate that the magnitude of the direct effect equals that of the indirect effect. We may formally run a test to verify it.
  • The LR test is \(\chi^2(df=1) = 0.004, p=.95\) suggesting that there is no evidence to reject the hypothesis of the direct effect equals the indirect effect.
## Proposed model
model <- "Gratitude ~ cp*Mindfulness + b*EI
          EI ~ a*Mindfulness
          Mindfulness ~~ 1*Mindfulness
          ## Define direct, indirect, and total effects
          Direct := cp
          Indirect := a*b
          Total := a*b + cp
          ## Apply a constraint on the direct and indirect effects
          cp == a*b"

RAM2 <- lavaan2RAM(model, obs.variables = varlist)
RAM2
## $A
##             Mindfulness EI      Gratitude
## Mindfulness "0"         "0"     "0"      
## EI          "0.1*a"     "0"     "0"      
## Gratitude   "0.1*cp"    "0.1*b" "0"      
## 
## $S
##             Mindfulness EI             Gratitude                   
## Mindfulness "1"         "0"            "0"                         
## EI          "0"         "0.5*EIWITHEI" "0"                         
## Gratitude   "0"         "0"            "0.5*GratitudeWITHGratitude"
## 
## $F
##             Mindfulness EI Gratitude
## Mindfulness           1  0         0
## EI                    0  1         0
## Gratitude             0  0         1
## 
## $M
##   Mindfulness EI Gratitude
## 1           0  0         0
## 
## $mxalgebras
## $mxalgebras$Direct
## mxAlgebra 'Direct' 
## $formula:  cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Indirect
## mxAlgebra 'Indirect' 
## $formula:  a * b 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$Total
## mxAlgebra 'Total' 
## $formula:  a * b + cp 
## $result: (not yet computed) <0 x 0 matrix>
## dimnames: NULL
## 
## $mxalgebras$constraint1
## MxConstraint 'constraint1' 
## $formula:  cp == a * b
stage2b <- tssem2(stage1, RAM=RAM2, intervals.type = "LB")
summary(stage2b)
## 
## Call:
## wls(Cov = pooledS, aCov = aCov, n = tssem1.obj$total.n, RAM = RAM, 
##     Amatrix = Amatrix, Smatrix = Smatrix, Fmatrix = Fmatrix, 
##     diag.constraints = diag.constraints, cor.analysis = cor.analysis, 
##     intervals.type = intervals.type, mx.algebras = mx.algebras, 
##     mxModel.Args = mxModel.Args, subset.variables = subset.variables, 
##     model.name = model.name, suppressWarnings = suppressWarnings, 
##     silent = silent, run = run)
## 
## 95% confidence intervals: Likelihood-based statistic
## Coefficients:
##    Estimate Std.Error   lbound   ubound z value Pr(>|z|)
## a  0.401476        NA 0.338588 0.465932      NA       NA
## b  0.270746        NA 0.215989 0.330137      NA       NA
## cp 0.108698        NA 0.087826 0.130046      NA       NA
## 
## mxAlgebras objects (and their 95% likelihood-based CIs):
##                   lbound Estimate    ubound
## Direct[1,1]   0.08782637 0.108698 0.1300457
## Indirect[1,1] 0.08782431 0.108698 0.1299560
## Total[1,1]    0.17576707 0.217396 0.2599277
## 
## Goodness-of-fit indices:
##                                                 Value
## Sample size                                13497.0000
## Chi-square of target model                     0.0052
## DF of target model                             0.0000
## p value of target model                        0.0000
## Number of constraints imposed on "Smatrix"     0.0000
## DF manually adjusted                           0.0000
## Chi-square of independence model             249.6905
## DF of independence model                       3.0000
## RMSEA                                          0.0000
## RMSEA lower 95% CI                             0.0000
## RMSEA upper 95% CI                             0.0000
## SRMR                                           0.0015
## TLI                                              -Inf
## CFI                                            1.0000
## AIC                                            0.0052
## BIC                                            0.0052
## OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
## Other values indicate problems.)
## Compare the models with and without c=a*b
anova(stage2, stage2b)
##                 base         comparison ep     minus2LL df      AIC      diffLL
## 1 TSSEM2 Correlation               <NA>  3 2.893374e-24 -3 6.000000          NA
## 2 TSSEM2 Correlation TSSEM2 Correlation  3 5.222545e-03 -2 6.005223 0.005222545
##   diffdf         p
## 1     NA        NA
## 2      1 0.9423893

Tips in MASEM

It takes a long time to run large models!

  • We may speed up the analyses by using more cores available in your computers.
library(metaSEM)

## Use all cores-2, i.e., keep 2 cores for other things.
mxOption(key='Number of Threads', value=parallel::detectCores()-2)

How do I know if my correlation matrices are “proper”?

  • There are restrictions on correlation matrices, e.g., they cannot larger than 1, which is known as positive definite. The results may be problematic when some of the correlation matrices are non-positive definite.
  • We may check the definite of the matrices with is.pd() function.
## X1 is okay
X1 <- matrix(c(1, 0.8, 0.8, 1), nrow=2)
X1
##      [,1] [,2]
## [1,]  1.0  0.8
## [2,]  0.8  1.0
is.pd(X1)
## [1] TRUE
## X2 is not okay as the correlation is 1.2.
X2 <- matrix(c(1, 1.2, 1.2, 1), nrow=2)
X2
##      [,1] [,2]
## [1,]  1.0  1.2
## [2,]  1.2  1.0
is.pd(X2)
## [1] FALSE
## X3 is not okay as the correlation matrix is not symmetric.
X3 <- matrix(c(1, 0.6, 0.8, 1), nrow=2)
X3
##      [,1] [,2]
## [1,]  1.0  0.8
## [2,]  0.6  1.0
isSymmetric(X3)
## [1] FALSE

What should I do when the solution is non-convergent?

  • Statistical packages use iterative procedures to estimate the parameters. When the iterations are convergent (code 0 or 1 in OpenMx), we can trust the results; otherwise, we should not interpret the results.
  • The results are okay:
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values indicate problems.)
  • The results are not okay:
OpenMx status1: 5 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)
Warning message:
In print.summary.meta(x) :
  OpenMx status1 is neither 0 or 1. You are advised to 'rerun' it again.

## Suppose the results are stored in an object called "fit"
fit <- rerun(fit)
summary(fit)
  • A common reason why there are non-convergent solutions is that there is not enough data. We may use the following code to check it.
pattern.na(Hunter83$data, show.na=FALSE)
##               Ability Job_knowledge Work_sample Supervisor
## Ability            13            11          11         12
## Job_knowledge      11            12          10         11
## Work_sample        11            10          12         11
## Supervisor         12            11          11         13

Exercises

Further resources

Thank you!

  • You have survived this workshop!
https://m.xkcd.com/1447/
https://m.xkcd.com/1447/

  1. Browne, M. W., & Cudeck, R. (1993). Alternative ways of assessing model fit. In K. A. Bollen & J. S. Long (Eds.), Testing Structural Equation Models (pp. 136-162). Newbury Park, CA: Sage.↩︎

  2. Hu, L., & Bentler, P. M. (1999). Cutoff criteria for fit indices in covariance structure analysis: Conventional criteria versus new alternatives. Structural Equation Modeling, 6, 1-55.↩︎

  3. Marsh, H. W., & Hau, K. T., & Wen, Z. (2004). In search of golden rules: Comment on hypothesis-testing approaches to setting cutoff values for fit indices and dangers in overgeneralizing Hu and Bentler’s (1999) findings. Structural Equation Modeling, 11, 320-342.↩︎

  4. Guyatt, G. H., Sackett, D. L., Sinclair, J. C., Hayward, R., Cook, D. J., Cook, R. J., Bass, E., Gerstein, H., Haynes, B., Holbrook, A., Jaeschke, R., Laupacls, A., Moyer, V., & Wilson, M. (1995). Users’ Guides to the Medical Literature: IX. A Method for Grading Health Care Recommendations. JAMA, 274(22), 1800–1804. https://doi.org/10.1001/jama.1995.03530220066035↩︎

  5. Schmidt, F. L., & Hunter, J. E. (2015). Methods of meta-analysis: Correcting error and bias in research findings (3rd ed.). Sage.↩︎

  6. Hedges, L. V., & Olkin, I. (1985). Statistical methods for meta-analysis. Orlando, FL: Academic Press.↩︎

  7. Hedges, L. V., & Vevea, J. L. (1998). Fixed- and random-effects models in meta-analysis. Psychological Methods, 3(4), 486–504. https://doi.org/10.1037/1082-989X.3.4.486↩︎

  8. Hedges, L. V., & Vevea, J. L. (1998). Fixed- and random-effects models in meta-analysis. Psychological Methods, 3, 486-504.↩︎

  9. Borenstein, M., Hedges, L. V., Higgins, J. P. T., & Rothstein, H. R. (2009). Introduction to meta-analysis. John Wiley & Sons.↩︎

  10. Higgins, J. P. T., & Thompson, S. G. (2002). Quantifying heterogeneity in a meta-analysis. Statistics in Medicine, 21, 1539-1558.↩︎

  11. Borenstein, M., Higgins, J. P. T., Hedges, L. V., & Rothstein, H. R. (2017). Basics of meta-analysis: I2 is not an absolute measure of heterogeneity. Research Synthesis Methods, 8(1), 5–18. https://doi.org/10.1002/jrsm.1230↩︎

  12. Rubio-Aparicio, M., López-López, J. A., Viechtbauer, W., Marín-Martínez, F., Botella, J., & Sánchez-Meca, J. (2020). Testing categorical moderators in mixed-effects meta-analysis in the presence of heteroscedasticity. The Journal of Experimental Education, 88(2), 288–310. https://doi.org/10.1080/00220973.2018.1561404↩︎

  13. Konstantopoulos, S., & Hedges, L. V. (2004). Meta-analysis. In D. Kaplan (Ed.), Handbook of quantitative methodology for the social sciences (pp. 281-297). New York: Sage.↩︎

  14. Cheung, M. W.-L. (2021). Meta-analytic structural equation modeling. In Oxford Research Encyclopedia of Business and Management. Oxford University Press. https://doi.org/10.1093/acrefore/9780190224851.013.225↩︎

  15. National Research Council (1992). Combining information: Statistical issues and opportunities for research. Washington, D.C.: National Academy Press.↩︎

  16. Brown, S. P., & Stayman, D. M. (1992). Antecedents and consequences of attitude toward the ad: A meta-analysis. Journal of Consumer Research, 19, 34-51.↩︎

  17. Premack, S. L., & Hunter, J. E. (1988). Individual unionization decisions. Psychological Bulletin, 103, 223-234.↩︎

  18. Norton, S., Cosco, T., Doyle, F., Done, J., & Sacker, A. (2013). The Hospital Anxiety and Depression Scale: A meta confirmatory factor analysis. Journal of Psychosomatic Research, 74(1), 74-81.↩︎

  19. Viswesvaran, C., & Ones, D. S. (1995). Theory testing: Combining psychometric meta-analysis and structural equations modeling. Personnel Psychology, 48(4), 865-885.↩︎

  20. Becker, B. J. (1992). Using results from replicated studies to estimate linear models. Journal of Educational Statistics, 17(4), 341-362.↩︎

  21. Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29-40. http://doi.org/10.3758/s13428-013-0361-y↩︎

  22. Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40-64. http://doi.org/10.1037/1082-989X.10.1.40↩︎

  23. Jak, S., & Cheung, M. W.-L. (2020). Meta-analytic structural equation modeling with moderating effects on SEM parameters. Psychological Methods, 25(4), 430–455. https://doi.org/10.1037/met0000245↩︎

  24. Cudeck, R. (1989). Analysis of correlation matrices using covariance structure models. Psychological Bulletin, 105(2), 317-327.↩︎

  25. Cudeck, R. (1989). Analysis of correlation matrices using covariance structure models. Psychological Bulletin, 105(2), 317–327. https://doi.org/10.1037/0033-2909.105.2.317↩︎

  26. Jak, S., & Cheung, M. W.-L. (2020). Meta-analytic structural equation modeling with moderating effects on SEM parameters. Psychological Methods, 25(4), 430–455. https://doi.org/10.1037/met0000245↩︎

  27. Jak, S., & Cheung, M. W. (2022, March 16). Can findings from meta-analytic structural equation modeling in management and organizational psychology be trusted?. https://psyarxiv.com/b3qvn.↩︎

  28. Cheung, M. W.-L., & Chan, W. (2005). Meta-analytic structural equation modeling: A two-stage approach. Psychological Methods, 10(1), 40–64. https://doi.org/10.1037/1082-989X.10.1.40↩︎

  29. Cheung, M. W.-L. (2014). Fixed- and random-effects meta-analytic structural equation modeling: Examples and analyses in R. Behavior Research Methods, 46(1), 29–40. https://doi.org/10.3758/s13428-013-0361-y↩︎

  30. Digman, J. M. (1997). Higher-order factors of the Big Five. Journal of Personality and Social Psychology, 73(6), 1246-1256. http://doi.org/10.1037/0022-3514.73.6.1246↩︎

  31. Nohe, C., Meier, L. L., Sonntag, K., & Michel, A. (2015). The chicken or the egg? A meta-analysis of panel studies of the relationship between work–family conflict and strain. Journal of Applied Psychology, 100(2), 522–536. https://doi.org/10.1037/a0038012↩︎

  32. Schutte, N. S., Keng, S.-L., & Cheung, M. W.-L. (2021). Emotional intelligence (EI) mediates the connection between mindfulness and gratitude: A meta-analytic structural equation modeling study. Mindfulness, 12(11), 2613–2623. https://doi.org/10.1007/s12671-021-01725-2↩︎

  33. Cheung, M. W.-L. (2021). Meta-analytic structural equation modeling. In Oxford Research Encyclopedia of Business and Management. Oxford University Press. https://doi.org/10.1093/acrefore/9780190224851.013.225↩︎

  34. Cheung, M. W.-L. (2022). Synthesizing indirect effects in mediation models with meta-analytic methods. Alcohol and Alcoholism, 57(1), 5–15. https://doi.org/10.1093/alcalc/agab044↩︎

  35. Cheung, M. W.-L., & Hafdahl, A. R. (2016). Special issue on meta-analytic structural equation modeling: Introduction from the guest editors. Research Synthesis Methods, 7(2), 112–120. https://doi.org/10.1002/jrsm.1212↩︎

  36. Cheung, M. W.-L., & Hong, R. Y. (2017). Applications of meta-analytic structural equation modelling in health psychology: Examples, issues, and recommendations. Health Psychology Review, 11(3), 265–279. https://doi.org/10.1080/17437199.2017.1343678↩︎

  37. Valentine, J. C., Cheung, M. W.-L., Smith, E. J., Alexander, O., Hatton, J. M., Hong, R. Y., Huckaby, L. T., Patton, S. C., Pössel, P., & Seely, H. D. (2022). A primer on meta-analytic structural equation modeling: The case of depression. Prevention Science, 23(3), 346–365. https://doi.org/10.1007/s11121-021-01298-5↩︎